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I* 


Introduction: 

The  purpose  of  this  grant  was  to  develop  methods  to  estimate  the  mechanical  properties 
of  tissue,  especially  the  elasticity,  in  vivo.  These  methods  will  be  used  to  evaluate  those 
mechanical  properties  as  indicators  of  breast  cancer.  MR  images  were  used  to  measure  the 
displacement  of  tissue  resulting  from  a  low  frequency  vibration  of  the  tissue.  We  developed 
methods  of  vibrating  tissue  in  the  magnetic  field  of  the  MR  system,  developed  the  MR  pulse 
sequences  that  measure  the  resulting  displacement  and  developed  the  an  algorithm  to  reconstruct 
the  elasticity  from  the  measured  displacements.  The  algorithm  to  reconstruct  the  elasticity 
allows  resolutions  a  factor  of  at  least  ten  better  than  previous  methods. 
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Body: 

•  Technical  Objective  1:  Refine  the  MR  measurement  of  three  dimensional  displacement 
during  vibration.  Build  the  high  power,  low  noise  apparatus  to  measure  the  three 
dimensional  displacement  during  vibration.  Compare  the  measured  displacement  to 
measurements  from  a  calibrated  hydrophone  and  tune  the  system. 

We  have  developed  a  driver  to  vibrate  tissue  at  low  frequencies  in  the  MR  magnet.  It  is  made  of 
three  stacks  of  piezoelectric  crystals  that  produce  very  high  forces  and  have  a  very  linear  response. 
See  Fig.  1  below.  The  apparatus  is  described  in  the  manuscript  submitted  to  Medical  Physics  [2]  and 
included  in  the  Appendix.  We  have  measured  displacements  in  the  magnet  using  a  gradient  echo 
pulse  sequence  that  measures  displacements  in  all  three  directions.  The  measured  displacements 
agree  well  with  independent  measurements  using  a  Capacitec  Model  410-SC  probe.  The 
displacements  measured  with  the  MR  and  the  Capacitec  probe  agree  very  well;  the  gel  at  the  edges 
of  the  phantoms  are  within  25%  with  the  Capacitec  measurements  of  the  motion  of  the  plastic  wall 
holding  the  phantom. 


Figure  1:  The  measured  displacements  vs.  driving  voltage  at  several  frequencies. 
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We  also  developed  noise  reduction  methods  to  improve  the  MR  estimates  of  displacement 
[3,4,5,8,10]  especially  at  very  low  displacement. 

•  Technical  Objective  2:  Compare  the  elasticity  calculated  from  the  MR  displacements  with  a 
linear  model  to  known  elasticity’s  for  isotropic  and  anisotropic  materials.  Modify  the  scale  with 
the  vernier  to  measure  the  elasticity  manually  with  a  mechanical  method.  Build  phantoms  and 
compare  the  elasticity  calculated  from  the  MR  displacements  with  a  linear  model  to  known 
elasticity’s  for  isotropic  and  anisotropic  materials. 

We  have  implemented  finite  element  code  that  models  the  displacement  of  a  variety  of  gels  and 
tissue.  It  is  clear  from  the  MR  measurements  and  from  the  simulations  that  the  full  partial 
differential  equation  is  needed  to  accurately  describe  the  system.  Representative  examples  of  those 
simulations  are  shown  in  [1, 2,6,9].  Copies  of  the  posters  and  papers  are  in  the  Appendix.  We  have 
measured  the  elasticity  of  gel  phantoms  with  mechanical  methods,  by  estimating  the  wavelength  of 
100  Hz  vibrations  in  the  MR  displacement  images  and  by  the  inversion  of  the  partial  differential 
equations  using  the  method  we  developed  and  have  reported  on  in  Mag.  Reson.  Med.  [1].  The  three 
methods  agree  to  a  reasonable  extent.  However,  we  believe  the  best  way  to  estimate  the  elasticity  of 
homogeneous  gels  is  by  estimating  the  wavelength  in  the  MR  displacement  images.  Any  frequency 
effects  are  accounted  for  and  Poisson’s  ratio  need  not  be  estimated.  These  results  are  also  shown  in 
the  four  publications  listed  above  and  included  in  the  Appendix. 

•  Technical  Objective  3:  Establish  the  limits  of  the  linear,  lossless  elastic  model  of  tissue  motion 
during  vibration.  Find  the  dependence  of  MR  elasticity  on  the  frequency  and  amplitude  of  the 
vibration.  Measure  motion  perpendicular  to  the  direction  of  forced  vibration.  Measure  viscous 
losses  by  the  attenuation  of  the  vibration  across  the  material.  Develop  and  test  finite  element 
code  to  calculate  the  displacement.  Compare  the  measured  displacements  with  the 
displacements  calculated  with  a  finite  element  analysis  for:  Phantoms  with  simple  geometry. 
Complicated  phantoms.  Lean  meat  (probably  a  roast).  Meat  with  fat  and  muscle  (probably  a 
slab  of  bacon). 

We  have  established  the  limits  of  the  linear  model  with  simulations  and  measurements  on  phantoms 
(results  of  technical  objective  2)  without  getting  to  more  complicated  structures;  the  full  partial 
differential  equations  (PDE’s)  are  required  to  dscribe  motion  accurately.  We  have  measured  viscous 
losses  with  reconstruction  methods  [2,7,9].  We  have  developed  finite  element  code  to  calculate  the 
displacements  with  known  materials  (the  forward  problem)  and  compared  them  to  measured 
displacements  [1,2,6,7,9]. 

The  results  we  now  have  suggest  that  we  need  to  measure  steady  state  vibration  in  the  tissue  to 
allow  accurate  reconstruction  of  the  elasticity  using  the  PDE’s  for  harmonic  motion.  The  PDE’s  for 
harmonic  motion  assume  steady  state  motion.  It  is  possible  to  solve  the  PDE’s  for  transient  motion 
as  well  but  it  is  significantly  easier  to  solve  the  steady  state  equations  than  the  equations  for  transient 
motion.  We  have  designed  a  system  to  drive  the  tissue  in  steady  state  using  the  MR  system’s  clock 
so  the  vibration  remains  in  sync  with  the  MR  over  the  entire  acquisition  [11].  We  have  begun  to 
gather  steady  state  data  on  gels  and  the  reconstructions  are  very  good  [11]. 

We  have  made  some,  limited  tissue  measurements  as  well  as  measurements  on  gels.  We  were 
lucky  enough  to  get  some  breast  tissue  from  a  breast  reduction  surgery.  We  made  measurements  on 
that  tissue  instead  of  the  meat  we  proposed  because  it  represents  real  breast  tissue  much  better  than 
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meat  which  is  mostly  muscle.  Recent  publications  such  as  [Krouskop,  T.A,  et.  al.,  “Elastic  moduli 
of  breast  and  prostate  tissues  under  compression”  Ultrasonic  Imaging  20(4):260-74, 1998],  suggest 
that  there  are  significant  differences  between  breast  tissue  and  muscle.  The  MR  images  of  breast 
tissue  from  breast  reduction  surgery  are  excellent  but  the  motion  is  very  small  because  the 
mechanical  coupling  between  the  tissue  and  the  plexiglas  wall  is  poor  without  the  skin  to  hold  the 
tissue.  There  are  several  possible  approaches  to  the  coupling  problem  but  we  had  no  time  to  pursue 
them.  However,  we  did  make  a  first  effort  at  vibrating  breasts  of  volunteers.  We  have  measured 
adequate  vibrations  in  a  volunteer  but  the  apparatus  proved  to  be  too  uncomfortable  to  use.  We  are 
currently  modifying  the  apparatus  extensively  to  make  it  usable;  this  work  is  being  done  under  the 
new  NIH  grant. 

There  are  two  significant  aspects  to  our  reconstruction  algorithm  [1,2,7,11,12].  The  most 
significant  aspect  of  the  new  algorithm  is  that  it  reconstructs  the  elasticity  with  resolution  on  the 
order  of  the  pixel  size  of  the  MR  images.  Previous  elasticity  estimates  from  dynamic  displacements 
had  spatial  resolution  on  the  order  of  the  wavelength  of  the  vibration  which  is  much  larger  than  the 
pixel  size.  The  other  significant  aspect  of  the  algorithm  is  that  it  is  readily  extended  to  three 
dimensional  problems  [12].  The  zones  used  in  the  method  can  be  made  small  enough  that  they  can 
be  solved  in  three  dimensions.  There  are  several  groups  that  believe  the  three  dimensional  problem 
must  be  solved  to  get  accurate  estimates  of  the  elasticity  in  vivo. 

We  believe  that  extensive  in  vivo  measurements  should  be  possible  in  approximately  one  year. 
Clinical  results  should  then  yield  the  usefulness  of  elasticity  in  breast  cancer  detection.  The  current 
line  of  research  provides  a  workable  approach  to  generating  elasticity  from  dynamic  vibration.  The 
results  of  this  work  should  be  compared  to  the  elasticity  estimates  from  quasi-static  MR  and 
ultrasound  to  see  if  the  extra  information  generated  from  dynamic  measurements  is  productive. 
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Research  Accomplishments: 

•  Method  of  reconstructing  the  in  vivo  elasticity  from  measured  displacements.  The  method 
has  been  tested  in  2D  and  can  be  extended  to  3D. 

•  Developed  apparatus  to  vibrate  tissue  in  the  MR  with  large  forces  and  accurate 
synchronization  with  the  MR  imaging  system. 

•  New  noise  reduction  methods  were  developed  to  improve  the  MR  estimates  of  displacement. 


Reportable  Outcomes: 

•  Five  papers  and  five  peer  reviewed  presentations  with  abstracts  have  resulted  from  this  grant. 
The  papers  and  presentations  are  listed  in  the  Bibliography. 

•  The  methods  and  results  from  the  work  supported  by  this  grant  have  formed  one  section  of  a 
Program  Project  Grant  submitted  to  the  NIH  in  June  of  1998.  It  received  very  good  scores 
and  was  resubmitted  in  June  of  1999.  The  1999  resubmission  was  funded. 
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Conclusions: 

We  have  shown  that  partial  differential  equation  solutions  are  needed  to  accurately  estimate  the 
elasticity  from  measured  displacements,  even  with  very  small  displacements.  Further,  we  have 
proposed  and  tested  an  algorithm  to  solve  for  the  elasticity  from  measured  displacements.  The 
algorithm  performs  exceptionally  well  in  simulation  and  produces  generally  accurate  solutions 
from  measured  data.  We  have  good  evidence  that  most  of  the  inaccuracies  in  the  reconstruction 
of  measured  data  results  from  transient  effects  in  the  tissue  vibration.  The  most  significant 
aspect  of  the  new  algorithm  is  that  it  reconstructs  the  elasticity  with  resolution  on  the  order  of  the 
pixel  size  of  the  MR  images.  Previous  elasticity  estimates  from  dynamic  displacements  had 
spatial  resolution  on  the  order  of  the  wavelength  of  the  vibration  which  is  much  larger  than  the 
pixel  size. 

With  improved  steady  state  vibration,  which  we  have  developed  and  have  begun  testing,  and  the 
extension  of  our  reconstruction  algorithm  from  2D  to  3D,  which  we  have  started,  obtaining 
accurate  estimates  of  the  elasticity  of  tissue  in  vivo  is  possible.  Therefore,  clinical  evaluation  of 
elasticity  as  an  indicator  of  cancer  should  be  possible  within  a  year  or,  at  most,  two. 
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^  Magnetic  Resonance  in  Medicine  42:779-786  (1999) 

An  Overlapping  Subzone  Technique  for  MR-Based  Elastic 
Property  Reconstruction 

E.E.W.  Van  Houten,^*  K.D.  Paulsen,'' M.I.  Miga,''  F.E.  Kennedy,''  and  J.B.  'Weaver^ 


A  finite  element-based  nonlinear  Inversion  scheme  for  mag¬ 
netic  resonance  (MR)  elastography  is  detailed.  The  algorithm 
operates  on  small  overlapping  subzones  of  the  total  region  of 
interest,  processed  in  a  hierarchical  order  as  determined  by 
progressive  error  minimization.  This  zoned  approach  allows  for 
a  high  degree  of  spatial  discretization,  taking  advantage  of  the 
data-rich  environment  afforded  by  the  MR.  The  inversion  tech¬ 
nique  is  tested  in  simulation  under  high-noise  conditions  (15% 
random  noise  applied  to  the  displacement  data)  with  both 
complicated  user-defined  stiffness  distributions  and  realistic 
tissue  geometries  obtained  by  thresholding  MR  image  slices.  In 
both  cases  the  process  has  proved  successful  and  has  been 
capable  of  discerning  small  inclusions  near  4  mm  in  dia¬ 
meter.  Magn  Reson  Med  42:779-786, 1999.  ©  1999  Wiley-Liss, 
Inc. 

Key  words:  elasticity  reconstruction;  nonlinear  inversion;  finite 
element  method;  magnetic  resonance  elastography;  subzone 
technique;  model-based  imaging 

The  diagnostic  value  of  tissue  elasticity  has  long  been 
appreciated  in  a  broad  spectrum  of  medical  applications, 
and  understanding  of  its  importance  continues  to  grow. 
From  pathology  detection  (1,2),  to  robotic  surgery  (3,4),  to 
the  use  of  computational  modeling  during  surgical  proce¬ 
dures  (5—7),  a  demand  for  detailed  and  accurate  tissue 
elasticity  information  has  been  generated.  Recent  research 
into  the  use  of  mechanical  properties  of  biological  tissue 
for  clinical  decision  making  has  moved  away  from  direct 
mechanical  measurements  (8,9)  and  turned  toward  various 
medical  imaging  technologies  to  assess  tissue  behavior 
under  mechanical  loads.  The  idea  of  ultrasound  elastogra¬ 
phy  has  been  introduced  (10-16)  in  which  some  form  of 
ultrasonic  displacement  measurement  technique  is  used  to 
detect  subsurface  tissue  motion.  This  displacement  infor¬ 
mation  can  then  be  correlated  to  the  elastic  property 
distribution  in  the  tissue  with  the  aid  of  a  model  for  tissue 
motion  as  a  function  of  shear  or  Young’s  modulus  (17-19). 

Ultrasound’s  inherent  lack  of  lateral  resolution  and 
limited  axial  resolution  when  compared  with  other  clini¬ 
cally  available  imaging  modalities  has  motivated  the  devel¬ 
opment  of  magnetic  resonance  (MR)  elastography  methods 
(20-22).  MR  offers  the  potential  of  generating  highly 
resolved,  three-dimensional  (3D)  information  with  relative 
ease,  as  opposed  to  the  considerable  challenge  associated 
with  obtaining  equivalent  data  from  ultrasound  tech¬ 
niques.  However,  given  the  availability  of  finely  sampled 
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3D  displacement  fields,  the  task  of  developing  a  robust 
algorithm  capable  of  deducing  elastic  property  distribu¬ 
tions  from  these  displacement  images  while  maintaining 
the  refinement  of  the  MR  data  remains. 

Strategies  for  addressing  the  reconstruction  problem 
have  varied  widely  to  date,  Chenevert  et  al.  are  investigat¬ 
ing  MR  elastography  through  a  quasistatic  displacement 
approach  (23),  whereas  Raghavan  and  Yagle  have  devel¬ 
oped  an  inversion  technique  based  on  a  finite  difterence 
formulation  of  the  global  elasticity  equations  (24).  A  collabo¬ 
ration  between  Lewa  and  De  Certaines  has  attempted  to 
determine  the  viscoelastic  properties  directly  from  MR 
measurements  (25).  Manduca  et  al.  have  developed  an 
inversion  scheme  based  on  a  local  frequency  estimation 
that  is  correlated  to  a  local  elasticity  value  (26),  and  Sumi 
and  Nakayania  have  presented  a  method  for  numerically 
integrating  the  two-dimensional  (2D)  stress-strain  relations 
to  reconstruct  a  shear  modulus  distribution  from  strain 
measurements  (27), 

In  this  report,  we  present  a  finite  element— based  method 
for  solving  the  elastography  inversion  problem  by  use  of  a 
least-squares  minimization  of  the  difference  between  mea¬ 
sured  displacement  data  from  the  MR  and  computed 
displacement  solutions.  Our  approach  is  not  unlike  that 
recently  presented  by  Kallel  and  Bertrand  for  ultrasound 
techniques  (18),  except  that  model-based  optimization  is 
performed  on  small  overlapping  subzones  of  the  total 
tissue  region  of  interest  that  are  processed  in  an  hierarchi¬ 
cal  order  determined  by  progressive  error  minimization. 
This  is  a  significant  shift  in  the  conceptual  framework  for 
property  inversion  that  allows  the  recovery  of  an  elasticity 
distribution  at  the  MR  displacement  measurement  resolu¬ 
tion.  Property  estimation  at  the  MR  pixel  level  is  not 
computationally  viable  as  a  single  global  minimization 
problem;  however,  the  subzone  approach  we  have  identi¬ 
fied  eliminates  this  limitation  by  recasting  the  image 
reconstruction  objective  as  a  sum  of  minimizations  rather 
than  a  single  minimization  of  sums.  The  results  show  that 
the  overlapping  zone  concept  is  robust  with  respect  to 
simulated  measurement  noise  and  that  local  minimization 
of  the  least-squares  match  between  the  model  and  the  MR 
displacement  data  leads  to  high-quality  global  property 
distribution  images. 


SUBZONE  INVERSION 

Our  approach  capitalizes  on  recent  advances  in  model- 
based  image  reconstruction  of  tissue  properties  whereby 
the  nonlinear  relationship  between  the  physical  property 
distribution  to  be  determined  and  the  measured  tissue 
response  to  an  applied  stimulus  is  preserved  (28-30). 
Specifically,  we  formulate  the  MR  elastography  image- 
reconstruction  problem  as  a  constrained  optimization  task 
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FIG.  1 .  Schematic  diagram  of  the  subzone  concept.  Q,  total  problem 
domain;  T,  boundary;  single  subzone  domain;  single  subzone 
boundary. 


whose  objective  is  to  minimize  the  difference  between  a  set 
of  measured  displacement  fields  and  those  computed  by  a 
model  description  in  which  the  tissue  property  distribu¬ 
tion  is  parameterized  as  a  set  of  unknown  coefficients.  The 
typical  strategy  is  to  define  a  single  objective  to  be  mini¬ 
mized  that  is  the  sum  of  the  squared  differences  between 
measured  and  calculated  quantities  over  the  entire  set  of 
tissue  response  observations  that  are  available  (31): 

min  F(E),  [la] 


FziEz),  for  the  zth  subzone.  For  Q  subzones, 

F(E)  =  2  [2] 

Z=1 


where  the  minimization  of  the  sum  is  replaced  by  the  sum 
of  minimizations  on  the  individual  subzones: 


min  F[E)  -  min 


•  2  min  F^jlFJl  [3a] 


with 

F,(F,)  =  E  -  “fel-  + 

when  each  subzone  consists  of  Nz  tissue  response  observa¬ 
tions  and  Mz  tissue  property  parameters  such  that  Nz  N 
and  Mz  M. 

The  advantages  of  this  approach  are  several  fold.  First, 
the  nonlinear  minimization  process  occurs  on  only  Mz 
tissue  property  parameters  using  Nz  observations  at  a  time. 
The  significant  reduction  in  the  size  of  the  inversion 
problem  is  important  because  the  least-squares  approach 
scales  cubically  in  the  number  of  optimization  parameters 
to  be  determined.  Second,  it  maximizes  the  utilization  of 
the  complete  MR  displacement  data  set  and  the  concomi¬ 
tant  tissue  property  resolution  that  can  be  achieved.  Assum¬ 
ing  that  tissue  displacements  can  be  measured  at  the  MR 
pixel  level,  the  total  amount  of  tissue  response  data  and 
tissue  property  values  that  could  be  recovered  in  a  single 
minimization  problem  exceeds  the  computational  re¬ 
sources  available  today.  By  dividing  the  problem  into 
subzones,  high-resolution  (MR  pixel-level)  property  maps 
can  be  deduced  that  take  full  advantage  of  the  high  density 
of  tissue  measurements  that  the  MR  technique  provides. 

Once  defined  on  the  subzone,  the  minimization  problem 
proceeds  in  standard  fashion.  Determination  of  the  sub¬ 
zone  elastic  properties  requires  differentiation  of  Fz  in  Eq. 
[3b]  with  respect  to  each  of  the  Mz  property  parameters 
contained  in  Ez,  which  produces  the  nonlinear  system 


where 


F{E)  =  E  K  -  f'lV  +  (v}"  -  lib] 

/-I 

and  n}"  and  are  the  .y  and  y  vector  components  of  the 
measured  displacement  at  location  /,  while  and  are 
the  calculated  vector  components  at  the  same  position,  for 
a  total  of  N  different  locations.  E  is  the  M-dimensional 
vector  of  elasticity  parameters  that  is  expanded  on  a 
continuous  basis  sot,  cl>,  to  define  the  tissue  property 
distribution  of  the  total  region  of  interest,  il. 

This  total  problem  domain,  11,  may  be  thought  of  as  the 
union  of  multiple  “siihzones,”  il^,,  of  the  total  ROI,  as 
illustrated  in  Fig.  1,  so  that  we  may  rewrite  the  global 
functional.  F[E),  as  a  sum  of  locally  defined  functionals, 


'  f'F,.  4=1 

fUlf 


Yz  fUly 


=  0 


ilF,  ^  ''“4 
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Solution  of  this  equation  set  by  Newton’s  method  leads  to 
iterative  improvements  in  the  elastic  property  profile  such 
that 


£;(n+l)  =  £:n  4.  [53] 

where  is  the  property  update  vector  determined  from 
the  solution  of  the  regularized  matrix  system 


[(HJ  +  al)l  1A£J  =  |-/;1  .  I5bl 


on  the  same  finite  element  discretization: 


- |V  •  GVu  +  V(X  +  G)V  •  u  =  -  pwM 

bE, 


bG  bu  b 

=  V  ■  —  Vu  +  V  •  GV  —  +  V  —  (\  +  G)V  ■  u 

bE,.  bE,^  riE,. 

flu  flu 

+  V(\  +  G)V  .  —  =  -  pw==  —  .  [7a] 


which  when  rewritten  in  the  form 


with  /z  =  \f[\  fl . fhV'  and  the  approximate  Hessian 

matrix  (second-order  terms  are  dropped),  having  the 
elements 


</ 


DE^l 


dlf;  dlfj  dv^l  ()v^ 

‘x  ‘z  ^  ‘z  z 

dE^l  dEl  <)El  c)Ei 


[5c] 


evaluated  at  the  current  property  estimate  denoted  by 
the  superscript  n.  In  Eq.  (5b),  a  is  a  scalar  regulari¬ 
zation  parameter  added  to  the  diagonal  of  H  to  facilitate 
its  inversion,  because  H  is  known  to  be  poorly  condi¬ 
tioned.  This  parameter  is  scaled  to  the  subzone  prob¬ 
lem  at  hand  by  use  of  the  Levenberg-Marquardt  approach 
(32). 

Solution  of  Eq.  [5b]  requires  a  vehicle  for  calculating  the 
subzone  displacement  field  and  its  derivatives  with  re¬ 
spect  to  each  property  parameter  given  the  current  estimate 
of  the  property  distribution  on  the  subzone.  Here,  we 
assume  that  the  displacement  field  is  described  by  the 
partial  differential  equation  governing  time-harmonic,  iso¬ 
tropic,  linearly  elastic  motion: 


V  •  GVu  +  V(\  +  G)V  •  u  =  -  pw^u,  [6a] 

where  ii  is  the  displacement  vector,  p  is  the  tissue  density 
and  G  and  X  are  Lame’s  constants 


E 

^  ”  (2  +  2v) 

[6b] 

vE 

^  “  (1  +  v)(l  -  2vj 

[6c] 

for  Poisson’s  Ratio,  v,  and  Young’s  modulus,  E.  For  simplic¬ 
ity,  we  consider  p  and  v  to  be  known  constants,  leaving  G, 
or  equivalently  E,  as  the  elastic  property  parameter  distri¬ 
bution  to  be  estimated  from  the  displacement  field.  To 
solve  Eq.  [6a],  we  use  the  finite  element  method  as 
summarized  in  the  Appendix. 

The  required  derivatives  are  calculated  by  differentiat¬ 
ing  Eq.  [6a]  directly  with  respect  to  each  property  param¬ 
eter,  Ey  for  ;  =  1,2,,..,  My,  and  solving  the  resulting  partial 
differential  equation  in  the  derivative  quantity  of  interest 


fJG 

V .  GVu;  +  v(x  +  G)v .  u;  =  -  pw-u;  ^  ^ 

-  V  -^  (\  +  G)V  •  u,  [7b] 

()E^ 

where  u;  =  ciu/dEy.*  has  the  identical  form  of  Eq.  l6a]  in  the 
quantity  u;  except  for  the  occurrence  of  two  additional 
right-hand-side  quantities  expressed  in  terms  of  u.  Because 
Eq.  [5b]  is  evaluated  at  the  current  properly  estimate,  u  can 
be  computed  through  Eq.  [6al  leaving  uy  as  the  only 
unknown  in  Eq.  [7bl.  Evaluation  of  the  terms  expressing 
the  differentiation  of  the  elastic  property  distribution  with 
respect  to  its  parameterization,  dG/dEy.,  is  facilitated  by 
expanding  the  elastic  constants  in  the  finite  element  basis 
so  that  E  =  Ez4)/.  The  finite  element  discretization  of 

Eq.  [7b]  is  describecl  in  further  detail  in  the  Appendix. 

In  practice,  the  subzone  inversion  algorithm  begins  with 
an  initial  estimate  of  the  elastic  property  distribution, 
defined  over  the  entire  problem  space,  H.  From  this 
estimate,  a  global  displacement  field,  u^  is  computed  firom 
Eq.  [6a]  with  the  finite  element  method  based  on  known 
displacement  conditions  from  the  MR  data  set,  u"',  applied 
along  the  global  boundary,  F,  The  squared  error  between 
the  resulting  displacement  solution  and  the  measured  MR 
data  is  then  calculated  for  each  element.  By  using  this  error 
metric,  an  hierarchical  list  of  element  centroids  is  gener¬ 
ated  in  which  the  element  order  is  based  on  a  decreasing 
squared-error  contribution.  A  subzone  domain,  Hz,  is  then 
formed  about  an  element  centroid  in  the  list  by  including 
all  nearby  elements  whose  centroids  are  within  a  user- 
defined  radial  distance  from  the  subzone  center.  Figure  2 
shows  an  example  of  a  simple  global  finite  element  mesh 
with  a  close-up  view  of  a  single  subzone  that  has  been 
extracted  for  illustrative  purposes.  Once  the  subzone  has 
been  identified,  a  displacement  field  is  calculated  on  the 
subzone  by  using  the  latest  property  parameter  estimate, 
E"f  and  the  MR  displacement  information  on  the  subzone 
boundary,  Fz,  as  the  boundary  conditions  required  for 
finite  element  solution.  The  subzone  property  distribution 
is  iteratively  updated  with  the  inversion  process  embodied 
in  Eq.  [5b]  until  a  local  convergence  criterion  between  the 
computed  and  measured  displacements  internal  to  the 
subzone  has  been  reached.  At  this  point,  the  next  element 
centroid  in  the  error  contribution  list  having  participated 
in  the  fewest  inversion  operations  is  used  to  detine  another 
subzone,  and  the  process  of  local  convergence  in  the 
displacement  field  between  computed  and  measured  quan- 


Y  Position  (m] 


782 


Van  Houten  et  al. 


FIG.  2.  Illustration  of  the  mesh-based  subzone  technique  in  which  a 
single  subzone  (a)  has  been  extracted  from  a  global  finite  element 
mesh  (b). 


titles  is  repeated.  The  zoning  process  continues  until  every 
element  in  the  global  mesh  has  been  iterated  a  minimum 
number  of  times.  The  subzone  solutions  then  end,  another 
global  displacement  field  calculation  is  executed  with  the 
latest  property  profile,  and  the  zoning  procedure  begins 
again.  Figure  3  illustrates  the  overall  image  reconstruction 
process. 


RESULTS 

Initial  evaluations  of  the  zoning  algorithm  have  been 
performed  on  various  numerical  simulations  of  the  dis¬ 
placement  patterns  generated  in  vibrating  tissue.  For  these 
experiments,  synthetic  data  were  produced  with  the  finite 
element  method  displacement  calculation  described  in  the 
Appendix.  The  mesh  geometries  were  developed  from  MR 
images  taken  from  actual  patients,  one  being  a  modified 
anatomically  coronal  breast  slice  and  the  other  a  coronal 
brain  image.  For  simplicity,  we  have  assumed  that  the  onl\ 
unknown  elastic  parameter  is  Young’s  modulus,  although 
this  is  not  an  inherent  limitation  in  our  algorithmic  ap¬ 
proach  per  se.  Values  for  tissue  density  and  Poisson  s  ratio 
were  prescribed  as  1020  kg/cm^  and  0.49,  respectively.  The 
small  wavelengths  that  develop  in  soft  tissue  (a  Young’s 
modulus  of  8000  Pa  was  used  for  the  background  tissue 
value  here]  require  that  the  planar  MR  image  be  divided 
into  a  large  number  of  elements  to  ensure  that  the  wave 
propagation  is  adequately  well  resolved.  For  example,  the 
mesh  used  for  the  breast  geometry  consisted  of  16,555 
nodes  and  32,635  linear  triangular  elements,  resolving  the 
tissue  continuum  to  approximately  0.8  mm.  The  wave¬ 
length  of  a  100-Hz  shear  wave  in  this  case  is  1.62  cni/cycle. 
so  this  resolution  provides  approximately  18  or  19  nodes 
per  mechanical  wavelength.  For  synthetic  data  production, 
a  spatial  elasticity  distribution  is  required.  We  have  gener¬ 
ated  property  distributions  by  thresholding  MR  images  that 
contain  either  an  arbitrary,  user-defined  elasticity  map,  as 
shown  in  Figure  4a,  or  one  that  follows  tissue  substruc¬ 
tures  that  are  identifiable  in  the  MR  image,  as  represented 


FIG.  3.  Flowchart  of  the  subzone  inversion 
algorithm. 
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FIG.  4.  Breast  elasticity  probiem  consisting  of  an  artificial  property  distribution  containing  three  ^jctsfrand  3  sS 

backqround  a-  Exact  Young's  moduius  distribution  (kPa),  which  includes  three  heterogeneities,  designated  as  objects  1 , 2,  and  3.  b.  Synthe  c 
?d2Sn  dfsplcement  Lgnitude  [pm],  c:  Synthetic  y-direction  displacement  magnitude  [pm],  d:  Reconstructed  Youngs  modulus 
distribution  (kPa)  in  the  presence  of  15%  measurement  noise. 


in  Figure  5a.  Before  inversion,  a  certain  percentage  of 
random  noise  is  added  to  the  synthetic  data  to  simulate 
signal  degradation  that  will  occur  in  practice.  This  noise  is 
generated  by  scaling  the  average  displacement  vyith  a 
random  number  up  to  a  given  percentage  and  adding  or 
subtracting  that  value  from  the  original  displacement  at  a 
particular  node.  The  noisy  solution  is  then  used  by  the 
inversion  algorithm  described  above  as  the  measured  data 
set. 

For  the  breast  case,  an  elasticity  distribution  was  created 
that  provided  challenging  inclusion  geometries  as  well  as  a 
variety  of  inclusion  stiffnesses.  The  background  tissue 
stiffness  was  assigned  a  Young’s  modulus  of  8000  Pa, 
which  is  believed  to  be  at  the  lower  end  of  actual  soft  tissue 
values  (8,9,33).  Inclusion  stiffness  ranges  were  determined 
as  multiples  (2X  for  object  1,  5  x  for  object  2  and  10 X  for 
object  3  in  Fig.  4a)  of  the  background  stiffness  to  test  the 
contrast  resolution  of  the  numerical  algorithm.  The  stiffest 
of  the  inclusions  (object  3  in  Fig.  4a)  is  roughly  4  mm  in 
diameter,  representing  a  very  small  tumor  within  the 
tissue.  Once  the  stiffness  information  has  been  lormulated 
for  the  forward  problem,  boundary  conditions  of  100  Hz 


and  10  pni  sinusoidal  displacements  are  applied  and  the 
displacement  solution  is  generated,  shown  here  as  x  and  y 
displacement  magnitudes  in  Fig.  4b  and  c,  respectively. 
Note  the  complex  nature  of  this  displacement  pattern.  For 
the  small  wavelengths  expected  in  soft  tissue,  these  com¬ 
plex  displacement  fields  could  lead  to  difficulties  in 
generating  propertv  distributions  based  on  direct  interpre¬ 
tation  of  the  displacement  or  strain  image.  The  property 
inversion  shown  in  Fig.  4d  was  generated  by  using  syn¬ 
thetic  data  with  15%  added  noise  with  an  initial  guess  of  a 
uniform  Young  s  modulus  of  7000  Pa.  The  inversion  pro¬ 
cess  consisted  of  18  sweeps  over  the  entire  space,  each 
sweep  involving  roughly  1000  subzones  of  approximately 
150  elements  and  100  nodes  each  to  ensure  that  every  node 
was  operated  on  at  least  once. 

To  demonstrate  the  ability  of  the  algorithm  to  process 
realistic  geometries  generated  from  MR  images,  an  elastic¬ 
ity  distribution  was  reconstructed  on  a  mesh  generated 
from  a  coronal  brain  slice  with  two  simulated  inclusions  as 
shown  in  Fig.  5a.  To  ensure  accurate  resolution,  this  mesh 
consisted  of  19,446  nodes  and  38,332  linear  triangular 
elements.  For  this  tissue,  a  background  Young’s  modulus  of 
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FIG.  5.  Brain  elasticity  problem  consisting  of  a  property  distribution 
derived  from  white/gray  matter  segmentation  of  a  coronal  MR  image, 
a:  Exact  Young’s  modulus  distribution  (kPa)  of  the  white/gray  matter 
background  with  two  artificially  placed  stiff  anomalies  designated  as 
objects  1  and  2.  b:  Reconstructed  Young’s  modulus  distribution  (kPa) 
in  the  presence  of  15%  measurement  noise. 

8000  Pa  was  used  (33],  representing  the  white  matter, 
while  the  gray  matter  was  assumed  to  be  twice  as  stiff.  The 
larger  inclusion  (designated  as  object  1  in  Fig.  5a}  measures 
approximately  1  cm  and  was  assigned  a  stiffness  five  times 
that  of  the  white  matter,  whereas  the  smaller  inclusion 
(designated  as  object  2  in  Fig.  5a)  is  roughly  4  mm  in 
diameter  and  was  specified  as  being  an  order  of  magnitude 
stiffer  than  white  matter.  The  inversion  process  was  com¬ 
pleted  in  six  global  sweeps  of  1240  subzones  on  average, 
with  approximately  164  elements  and  100  nodes  in  each 
zone.  Figure  5b  shows  the  recovered  property  distribution 
in  the  presence  of  15%  added  noise;  which  compares  very 
favorably  with  the  exact  distribution  contained  in  Fig.  5a. 
A  small  degree  of  spatial  filtering  was  found  to  be  helpful 
in  achieving  a  convergent  solution  in  this  case  of  high 
noise  (34).  This  technique  works  to  average  the  local 
property  value  at  node  /  with  the  values  at  its  immediately 
adjacent  nodes  so  that  E‘r  = 

where  N,-  is  the  number  of  neighbor  nodes  connected  to 


node  i.  For  the  inversion  shown  in  Fig.  5b,  a  value  of  0.2 
was  used  for  0.  Note  that  no  spatial  filtering  was  used  (i.e., 
0  =  0)  during  the  inversion  shown  in  Fig.  4d. 


CONCLUSIONS 

The  image  reconstructions  presented  in  Figs.  4  and  5  show 
that  the  zoned  inversion  scheme  is  able  to  generate  accu¬ 
rate  Youngs  modulus  distribution  images  on  the  basis  of 
displacement  information  obtained  in  the  presence  of  high 
noise  (up  to  15%).  The  zoned  inversion  method  also  allows 
for  a  high  degree  of  parameter  discretization,  taking  full 
advantage  of  the  data-rich  nature  of  an  MR  displacement 
image.  This  high  resolution  allowed  the  inversion  routine 
to  discern  hard  inclusions  as  small  as  4  mm  in  diameter 
during  simulation.  Although  there  are  additional  complexi¬ 
ties  associated  with  the  inversion  of  actual  MR  data,  these 
promising  simulation  results  suggest  that  the  subzone 
technique  should  provide  a  powerful  framework  for  recov¬ 
ering  elasticity  distributions  with  MR  elastography.  Specifi¬ 
cally,  it  provides  a  computationally  viable  approach  that 
capitalizes  on  the  MR  measurement  density  to  yield  high- 
resolution  (pixel-level)  tissue  property  maps.  The  complex 
nature  of  the  displacement  fields  that  develop  in  soft  tissue 
makes  an  imaging  algorithm  that  exploits  modeling  con¬ 
cepts  essential.  Use  of  the  finite  element  method  allows  the 
incorporation  of  tissue  mechanics  into  the  analysis  of 
measured  displacement  patterns  so  that  the  complicated 
waveforms  inherent  in  multidimensional  elastic  systems 
may  be  taken  into  account. 

The  finite  element  inversion  process  also  provides  some 
important  avenues  for  dealing  with  noisy  data  and  more 
complicated  physical  models.  The  iterative  inversion 
scheme  detailed  here  is  amenable  to  a  total  variation 
minimization  process,  which  can  be  useful  in  reducing  the 
effects  of  noise  degradation  (35).  Work  is  also  under  way  to 
incorporate  a  Maxwellian  damping  term  into  the  inversion 
process  so  that  tissue-damping  effects  can  be  both  ac¬ 
counted  for  and  measured  in  the  same  manner  that  Young  s 
modulus  is  recovered  with  the  algorithm  detailed  here. 
Furthermore,  the  process  should  be  adaptable  to  modeling 
transient  motion  rather  than  using  a  steady-state  assump¬ 
tion  if  steady-state  data  sets  cannot  be  obtained  from  the 
MR.  In  summary,  this  zoned  finite  element  inversion 
technique  provides  a  powerful  method  for  deriving  highly 
resolved  Young’s  modulus  distribution  information  in  a 
way  that  is  robust  in  the  presence  of  high  noise  and 
adaptable  to  a  variety  of  modifications  that  could  improve 
its  performance  in  the  future. 


APPENDIX 

Here  we  describe  the  finite  element  formulation  of  the 
forward  and  inverse  problems  in  more  detail.  The  govern¬ 
ing  equation  of  linear  elasticity,  Eq.  [6al,  is  cast  in  terms  of 
the  Cauchy  stress  tensor,. ^Twhich  leads  more  readily  to  the 
2D  plane  stress  or  plane  strain  conditions  assumed  for  the 
forward  and  inverse  solutions  presented. 
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The  Forward  Problem 

If  steady-state  harmonic  motion  is  assumed,  the  real  valued 
displacement  vector  solution  can  be  represented  as  the  real 
part  of  a  complex-valued,  spatially  varying  displacement 
phasor  multiplied  by  the  complex  exponential 

u[x,  y,  t]  =  /?e|u(x,  [8] 

where  spatial  coordinates  x  and  y  and  complex  displace¬ 
ment  components  ii  and  v  define 


u(x.y)  =  |^|.  [91 

In  stress  tensor  notation,  the  harmonic  equilibrium  condi¬ 
tion  will  then  be  (36) 


To  facilitate  the  solution  of  this  set  of  equations 
for  complicated  geometries  containing  realistic  property 
distributions,  a  finite  element  discretization  method 
is  adopted  where  u,  the  approximate  solution  to  the 
displacement  magnitude  u,  is  expanded  on  a  set  of  locally 
active,  spatially  varying  Lagrangian  basis  functions,  4)/,  so 
that 

u  =  2  “y'f’r  '■  2)  '0-4>y’ 

hi  hi 

where  Uj  and  Vj  are  the  approximate  x  and  y  di¬ 
rected  displacements  at  the  /th  of  N  nodes.  Developing  a 
Galerkin  weak  form  of  Eq.  [10]  generates  the  system  of 
equations 

lAl|“l  =  lb!,  [13] 


V  • ./  =  -  poj“u  +  F,  [10] 

making  use  of  tlie  complex  harmonic  inertia,  =  -  pco- 

ue'“^  and  including  any  additional  body  forces  F  that  may 
be  present.  The  stress  tensor. /^can  be  defined  as 

.T'  =  iT^  +  jTy  + 

Tx  =  +  )  V 

Ty  =  iV-  +  JCTy  + 

Tx  =  ITxx  +  ]-r,y  +  tor,, 

and  for  2D  plane  strain  or  plane  stress  assumptions  the 
stress-strain  relations  for  elastic  solids  can  be  written  in  the 
matrix  form 


with  matrix  [A]  consisting  of  the  subelements 
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and  the  column  vectors  of  unknown  nodal  displacements 
and  right-hand-side  forcing  being  written  as 


^  =  |ul  =  [iq,  Vi,  U2.  V2 - -  KvP’. 

V 


X  •  [#./•  n4>/ds  -  (F4),)] 

for  j  and  j  running  from  1  to  N. 

The  Inverse  Problem 

For  the  inversion  problem,  the  basis  set  used  to  expand  the 
parameterized  elasticity  solution  is  taken  to  be  the  same 
finite  element  basis  set  used  to  expand  the  displacement 
solution  u.  The  calculation  of  the  du/dE  terms  needed  to 
generate  the  Hessian  matrix  and  the  right-hand-side  vector 
fin  Eq.  [5b]  is  then  achieved  directly  by  the  differentiation 
of  Eq,  [13],  the  discretized  version  of  Eq.  [10],  with  respect 


where  E  is  Young’s  modulus. 
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to  E,  leading  to 
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The  determination  of  the  elastic  property  distribution  in  heterogeneous  gel  samples  with  a  finite 
dement  iltsed  reconstruction  scheme  is  considered.  The  algorithm  operate.s  on  sma  1  overlappm 
subzones  of  the  total  field  to  allow  for  a  high  degree  of  spatial  discret.zat.on  while 
computational  tractability.  By  including  a  Ma.Kwellian-type  viscoelastic  property  in  tte 
ics  and  optimizins  the  spatial  distribution  of  this  property  in  the  same  mannei  as  elasticity,  . 
Yoiino’s  Lduluslmage  is  obtained  which  reasonably  reflects  the  true  distribution  within  the  gel. 
However  the  imase  lacks  the  clarity  and  accuracy  expected  based  on  simulation  experience.  P 
lini  investiaaVions  suggest  that  transient  effects  in  the  data  are  the  cause  of  a  sjgnifican 
mismatch  betwe'en  the  inversion  model,  which  assumes  steady-state  conditions,  and  the  actua 
displacements  as  measured  by  a  phase  contrast  MR  technique.  ©  2000  American  Association  f 
Physicists  in  M^^^//c7«^'.[S0094-2405 (00)01 901-5] 

Key  words;  MR  elastography.  experimental  elasticity  reconstruction,  viscoelastic  tissue  behavior, 
finite  element  inversion  techniques 


I.  INTRODUCTION 

Interest  in  accurate,  quantitative  evaluation  of  the  physical 
properties  of  biological  tissue  is  rapidly  increasing  as  the 
value  of  this  information  has  been  appreciated  in  such  appli¬ 
cations  as  lesion  detection,  medical  examination,  and  com¬ 
puter  assisted  surgical  procedures.  One  particular  area  of  in¬ 
terest  includes  the  assessment  of  the  elastic  properties  of 
tissue  through  the  analysis  of  noninvasive  image-based, 
strain  measurements.  Ultrasound  elastography  has  been  stud¬ 
ied  for  some  time.'"^  where  ultrasonic  imaging  is  used  to 
detect  subsurface  displacement  response  to  an  externally  ap¬ 
plied  displacement  source.  While  ultrasound  elastography 
continues  to  hold  promise  in  terms  of  providing  valuable 
tissue  data,  it  is  presently  limited  by  the  lack  of  lateral  reso¬ 
lution  in  ultrasonic  imaging  compared  to  other  available  mo¬ 
dalities.  This  potential  drawback  has  led  to  the  invesuption 
of  MR-based  strain  imaging  techniques.’  '  In  addition  to 
fully  three-dimensional  imaging  capabilities,  MR  offers 
higher  resolution  than  ultrasound,  making  it  a  very  appealing 
method  for  discerning  micron  level  displacement  fields  in 

biolosical  structures.  ,  .  .  ,, 

Re"«ardless  of  the  strain  imaging  technique,  it  is  generally 
aoreed  that  some  type  of  model-based  reconstruction  proce¬ 
dure  is  required  for  determination  of  elastic  property 
distributions  To  date  these  models  have  varied  widely  in 
XlexW:  ITom  simple  algorithms  relating  a  Local  Fre¬ 


quency  Estimation  (LFE)  to  stiffness  via  a  standard  wave 
equation  formulation'^  to  nonlinear  inversion  schemes  which 
operate  on  a  partial  differential  equation  (PDE).' '  "  While 

many  of  the  methods  proposed  to  date  have  shown  promise 
in  simulation,  there  has  been  limited  success  in  producing 
accurate  inversions  of  measured  dynamic  strain  data.  This 
points  to  a  variety  of  issues,  such  as  viscoelastic  behavior 
and  transient  effects,  which  need  investigation  in  order  to 
develop  successful,  robust  algorithms  for  elasticity  imaging 
through  interpretation  of  harmonic  displacement  fields. 

While  inversion  in  the  case  of  dynamic,  rather  than  quasi- 
static,  mechanical  motion  has  shown  the  need  for  a  high 
level  of  algorithmic  sophistication,  the  benefits  of  the 
dure  are  potentially  significant.  By  generating  strain  fiel  s 
through  dynamic  displacement  propagation,  elasticity  infor¬ 
mation  can  be  inferred  in  structures  inaccessible  by  means  of 
quasi-static  surface  actuation,  e.g.,  the  brain.  In  addition,  dy¬ 
namic  di.splacements  are  useful  for  estimating  tissue  prefer 
ties  other  than  elasticity,  such  as  dispersion  and  density.  This 
paper  examines  the  first  applications  of  a  subzone-based 
elastic  property  reconstruction  scheme,  with  an  additional 
Maxwellian-type  viscoelastic  term,  to  experimental  MR  data. 
This  algorithm  is  based  on  a  well-established  nonlinear  in¬ 
version  scheme.'-'*-’’  where  the  global  reconstruction  tech¬ 
nique  is  instead  performed  on  small,  local  “subzones  of  the 
tissue  region  of  interest  (ROI).  The  improvements  in  algo- 
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rithmic  performance  found  through  the  addition  of  viscoelas- 
tic  behavior  are  encouraging  but  not  complete.  Evidence  pre¬ 
sented  here  suggests  that  the  method’s  assumption  of  steady- 
state  motion  in  the  MR  generated  data  may  be  the  primary 
discrepancy  in  the  data-model  match. 

II.  INVERSION  METHOD 

The  nonlinear  elastic  property  reconstruction  scheme  that 
forms  the  core  of  the  subzone  inversion  technique  is  centere^d 
around  a  squared  error  minimization  formulation," 
where  the  agreement  between  the  calculated  displacements 
derived  froni  the  most  recent  property  estimate,  ii‘;  at  loca¬ 
tion  /,  and  the  measured  displacement  data  at  that  location, 
u'/' ,  is  sensed  by  the  functional 

F=i:  (u'/'-u;-)-. 

/=! 

for  N  measurement  locations.  Minimization  of  this  tunctional 
is  accomplished  by  equating  its  derivatives  with  respect  m 
the  reconstruction  parameters  to  zero,  generating  the  matrix 
system 

[(H+«I)]{Ar}-{f}=0,  (2) 

with  AT  representing  the  vector  of  parameter  updates,  f  be¬ 
ing  the  derivatives  of  F  in  Eq.  (1)  with  respect  to  the  param¬ 
eter  field  r,  and  H  being  the  approximate  Hessian  matrix. 
The  regularization  parameter  a  is  added  to  the  diagonal  of  H 
and  scaled  to  the  inversion  problem  at  hand  by  the 
Levenberg-Marquardt  method"^  to  facilitate  solution  of  Eq. 
(2). 

Two  elastodynamic  models  have  been  implemented  to 
date  in  order  to  provide  displacement  calculations  based  on 
the  estimated  parameter  distributions.  The  simplest  involves 
an  undamped  linear  elastic  relationship,  while  the  other  adds 
a  Maxwellian  viscoelastic  coefficient  to  account  tor  the  at¬ 
tenuation  of  displacement  in  tissue.  Expressed  in  terms  of  a 
'  ’  — - written  as" 


PDE  in  displacement,  u,  these  models  can  be 


V-CVu-hV(X-4-G)V-u  =  p-^.  w 

for  a  linear  elastic  tissue  of  density  p  having  properties  G  and 
\  and 

du  (?‘U  ,  . 

V-GVu+V(\-l-G)V-ii=f  — +  p-^.  W 

for  a  Maxwellian  viscoelastic  material  with  damping  coeffi- 
cient 

In  general,  the  inversion  problem  defined  in  Eq.  (2)  is 
intractable  for  highly  resolved  material  property  distributions 
as  computine  a  single  solution  update  scales  cubically  wit 
the  size  of  the  vector  F.  To  overcome  this  resolution  limita¬ 
tion,  a  subzone  inversion  scheme  has  been  developed 
which  recasts  the  global  enor  functional  in  Eq.  (1) 
of  functionals  defined  on  small  subzones  of  the  total  ROl, 

F.  <» 


^lO.  1.  Piezodectric  dispUirenicm  actuator  with  ittiage  A  showing  the  entire 
levice  includina  the  phantom  bos  on  tlie  left  and  image  B  showing  a  close 
ip  of  the  nine  parallel  driven  piezocrystal  stacks  attached  to  one  wall  ot  the 
•inid  phantom  container.  The  imparted  motion  is  parallel  to  the  long  axis  a 
the  piezoelectric  actuator  which  is  denoted  as  the  .v  direction  m  subsequent 


where  min(F)  is  replaced  by 

X  min(F.) 

:=l 

for  Q  subzones.  Further  clelails  on  the  mathematical  and 
computational  underpinnings  of  this  technique  are  described 
in  Ref  tq-  the  aluorithm  essentially  relying  on  foundations 
previouslv  developed  in  Refs.  18.  19.  In  numerical  simul.v 
tions  the  formulation  has  been  found  to  be  very  effective  in 
reconstructing  high  resolution  (MR  pixel  level)  parameter 
distributions  even  in  the  presence  of  high  noi.se  (up  to  15%  in 
the  undamped  linear  elastic  case),  as  reported  in  Ref.  -0.  The 
subzone-based  inversion  method  has  also  been  found  to  ex¬ 
ecute  relatively  quickly,  with  global  property  distribution  so¬ 
lutions  for  meshes  of  approximately  20000  nodes  being 
completed  in  a  few  hours.  Direct  solution  of  the  full  inverse 
problem  for  a  mesh  of  this  size  is  computationally  mfeasib  e 
on  workstations  currently  available.  Here,  we  describe  om 
initial  experience  in  applying  the  algorithm  to  dynamic  MR 
displacement  data  measured  in  gel  phantoms. 

III.  MR  DISPLACEMENT  IMAGING  TECHNIQUE 

To  obtain  displacement  measurements,  we  have  devel¬ 
oped  an  MR-compatiblc  mechanical  driver  which  can  vibrato 
tissue-like  gel  phantoms  in  a  strong  magnetic  field,  igufv 
shows  the  device  which  is  currently  in  use.  The  design  i^ 
centered  around  an  actuator  constructed  of  nine  piezoelectric 
stacks  wired  in  parallel.  Three  layers,  each  separate  >  sta 
bilizing  ceramic  plates  and  consisting  ol  three  piezostac 
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Motion 


Motion  Sensibring  GadcntJ 


Fio  '»  MR  motion  sensitizing  gradient  sequence  showing  that  the  duration 
of  mechanical  motion  is  directly  linked  to  the  data  acquisition  time  associ¬ 
ated  with  one  phase  encoding  event. 

\  ere  used  to  provide  adequate  displacement  under  dynarnic 
load.  Individually,  each  stack  is  capable  of  ±7  fxm  of  dis¬ 
placement  when  driven  by  a  lOO-Vpp  signal,  so  that  the  as¬ 
sembly  as  a  whole  can  ideally  produce  ±21  fim  of  motion. 
The  stacks  are  powered  by  a  wide  band  amplifier  which 
boosts  the  sinusoidal  driving  signal  generated  by  the  MR 
pulse  program.  Using  the  MR  pulse  sequence  to  generate  the 
drivin<’  signal  for  mechanical  vibration  allows  almost  any 
waveform°to  be  used  and  also  allows  highly  accurate  control 
:  f  the  phase  between  the  mechanical  vibration  and  the  mo¬ 
tion  encoding  gradient.  However,  this  also  limits  the  time  of 
mechanical  stimulation  to  the  signal  duration  of  the  motion 
sensitizing  gradients  as  described  below  (see  Fig.  2).  Inde¬ 
pendent  measurements  of  actuator  motion  using  a  capaci- 
tively  coupled  displacement  probe  have  demonstrated  that 
the  response  of  the  actuators  is  linear  in  the  applied  voltage 
for  a  given  frequency,  yielding  10-15  /xm  of  displacement  at 
100  Hz,  60  Vpp. 

The  motion  encoding  MR  gradient  sequence,  based  on  a 
■.imple  gradient  echo  sequence,  is  illustrated  m  Fig.  2.  The 
ima-ino  process  involves  three  cycles  of  mechanical  motion 
driven  by  signal  from  the  RF,  the  last  two  of  which  coincide 
with  the  sensitizing  gradients.  The  phase  in  the  MR  image  is 
recorded  for  different  phase  relationships  between  the  ap¬ 
plied  motion  and  the  encoding  gradients  to  calculate  the  har¬ 
monic  displacement  in  each  pixel  of  an  image.  Four  such 
phase  relationships  are  used  to  ensure  internal  consistency 
between  the  measured  phase  and  the  fit  to  the  expression  C 
-t-M  cos  <p.  where  the  amplitude,  M,  and  the  relative  phase, 
(p,  completely  characterize  the  harmonic  motion.  The  con¬ 
stant  phase  offset,  C,  has  no  bearing  on  the  measured  motion. 

■  MR  timing  and  resolution  parameters  for  these  image  acqui¬ 
sitions  are;  rR  =  300m.s;  T£:=58ms;  FOV=8cm;  256 
X  128  pixels;  10-mm  slice  thickness. 


IV.  EXPERIMENTAL  RESULTS 

Initially,  homogeneous  gels  were  used  to  validate  the  MR 
measurement  system.  Figure  3  shows  a  typical  example  of 
the  displacements  obtained  in  a  homogeneous  gel  phantom 
(5  cmX5cmX4cm),  where  the  100-Hz,  x-directed  actuation 
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was  applied  through  the  bottom  and  two  side  walls  of  the 
rigid  gel  container  (see  Fig.  1).  The  displacements  were  re¬ 
corded  in  both  the  horizontal  (.v)  and  vertical  (y)  directions 
such  that  they  form  a  two  component  vector  field  in  the 
imaging  plane  parallel  to  the  axis  of  mechanical  translation. 
Beewse  of  the  symmetry  in  this  case,  only  a  single  central 
plane  is  illustrated.  Distortional  (shear)  wave  fronts  are  evi¬ 
dent  with  a  wavelength  of  approximately  1  cm.  This  wave¬ 
length  is  consistent  with  the  known  Young's  modulus  of  the 
kPa,  as  measured  in  separate  mechanical  tests). 
Peak  displacements  measured  in  the  gel  at  the  walls  of  the 
box  were  approximately  10  ixm.  which  mirrors  the  data  col¬ 
lected  with  an  independent  measurement  probe. 

In  Fia.  3,  the  images  in  the  first  column  represent  data 
collected  for  .v-directed  (horizontal)  motion  while  the  second 
column  shows  the  corresponding  data  for  y-directed  (verti¬ 
cal)  motion.  The  first  row.  A,  shows  a  standard  giay.sealt  MR 
image  of  the  homogeneous  gel  phantom,  while  rows  B  and  C 
show  the  reconstructed  motion  parameters  M  and  <f.  respec¬ 
tively.  The  bottom  row,  D,  shows  the  overall  displacement 
field,  given  by  M  cos  <p  for  each  vector  displacement  compo¬ 
nent.’  It  is  interesting  to  note  that  there  is  a  small  discrepancy 
between  these  displacement  images  and  the  know'n  boundary 
conditions,  i.e.,  a  nonzero  displacement  along  the  bottom  of 
the  gel  phantom.  This  is  due  to  an  imperfection  in  an  early 
prototype  of  the  mechanical  driver,  where  the  horizontal 
plane  of  motion  did  not  exactly  correspond  to  the  horizontal 
plane  of  the  MR.  leading  to  a  uniform  y-directed  displace¬ 
ment  offset  within  the  gel.  . 

The  development  of  wavelike  behavior  is  readily  visible 
in  the  magnitude  (row  B)  and  displacement  fields  (row  D)  as 
expected.  However,  the  mechanics  of  the  wave  propagation 
are  fairly  complex,  even  for  a  homogeneous  gel  phantom  as 
illustrated  in  Fig.  3,  due  to  the  asymmetric  (top  to  bottom) 
boundary  conditions  which  exist.  From  the  .v  directional  dis¬ 
placement  (row  D,  left),  it  is  evident  that  the  predominant 
mode  of  wave  propagation  is  shear  motion  caused  by  the 
bottom  surface  acting  as  a  distortional  wave  source.  The  uni¬ 
form  horizontal  peaks  and  valleys  suggest  the  gel  is  nearly 
incompressible  at  this  length  scale  (short  compared  to  the 
dilatational  wavelength)  and  low  frequency,  so  that  dilata- 
tional  wave  propagation  is  not  supported.  The  same  conclu¬ 
sion  is  reached  by  examining  the  y  motion  in  Fig.  3  (row  D, 
right)  where  uniform  vertical  contours  exist,  indicating  a 
shear  mode  of  propagation.  Another  noticeable  characteristic 
is  the  apparent  damping  of  motion  from  the  walls  to  the 
center  of  the  gel  (row  B)  suggesting  that  the  wavefronts  may 
have  undergone  some  dispersion  or  may  not  represent  fully 
developed  harmonic  motion. 

Figure  4  shows  MR  generated  .v  and  y  displacement  data 
for  a  block  of  agar  gel  with  a  hard,  spherical  agar  gel  inclu¬ 
sion  approximately  1.5  cm  in  diameter.  The  plane  shown  is 
the  same  as  Fig.  3,  cutting  through  the  center  ol  the  phantom 
with  its  .v-axis  in  line  with  the  direction  of  vibration^  Ihe 
v-directed  displacement  offset  associated  with  Fig.  -  as 
been  removed  by  the  elimination  of  the  imperfection  in  the 
mechanical  driver.  Figure  5  presents  a  standaid  giay  sca  e 
MR  image  of  the  gel  phantom.  Interestingly,  the  inclusion  is 
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Fig.  3.  Raw  data  from  MR  motion  en¬ 
coding  gradients  using  a  honiogeneou> 
gel  phantom  within  the  device  shown 
\n  Fig.  I  with  voxel  intensity  grayscale 
images  (row  A),  magnitudes  of  .v  and  y 
displacement  respectively  [fxm\  (row 
B),  phase  of  A’  and  y  directional  dis¬ 
placement  [rad]  (row  C).  and  recon¬ 
structed  A  and  y  displacement  tro;*. 
magnitude  and  phase  [^m]  (row  D 
The  left-hand  column  represents  a  di¬ 
rected  motion  sen.sitivity  while  the 
right  represents  y  directed  sensitivity. 
Note  that  the  images  in  this  figure  are 
dimensioned  in  pixel  coordinates. 


B  10  12  -20 


Microns 


biuely  visible  in  image  A.  showing  that  the  elastic  property 
contrast  does  not  translate  into  conventional  MR  image  con¬ 
trast  It  is  also  the  case  that  the  displacement  fields  them¬ 
selves  (Fig.  4)  do  not  clearly  indicate  the  presence  ot  the 
harder  incKision  embedded  within  the  gel.  For  visualization 
purposes,  the  area  of  heterogeneity  has  been  outlined  in  im- 

a^e  B  ol  Fig.  5.  . 

"  Image  rcconsiruciions  of  the  elastic  property  distribution 

in  the  heteroiieneous  gel  sample  (Fig.  5)  obtained  trom  the 
MR  displacement  data  (Fig.  4)  are  presented  in  Fig.  6.  Two 
inversions  are  shown;  one  with  the  undamped  elastic  mode 
from  Eu.  (.1)  (image  A)  and  the  other  with  the  damped  mode 
presented  in  Eq.  (4)  (images  B  and  C).  where  image  B  shovvs 
the  elasticity  distribution  and  image  C  shows  the  map  of  the 
Maxwellian  alienuaiion  coeflicient  for  the  damped  model  re¬ 
construction.  While  both  images  A  and  B  clearly  recover  the 
presence  of  a  hard  inclusion  near  the  center  ol  the  gel.  the 
size  and  overall  localization  of  the  inclusion  are  more  accu- 
ratelv  represented  in  the  damped  model  inversion.  Image  C 
indicates  that  the  surrounding  gel  has  a  relatively  low  attciiu 


ation  component  while  the  inclusion  exhibits  a  level  of  in- 
creased  damping. 

The  results  in  Fig.  6  are  encouraging,  especially  since 
they  represent  one  of  the  first  elastic  property  images  dern  ed 
from  experimental  MR  data  where  displacements  have  been 
measured  under  conditions  of  a  dynamically  applied  me¬ 
chanical  stimulus.  However,  further  improvements  are  desir¬ 
able  and  necessary  if  MR  elastography  is  to  progress  clini¬ 
cally.  While  qualitatively  satisfying,  the  images  in  Fig.  6  are 
not  immediately  as  promising  as  those  achieved  in  simu  a- 
lion  studies,-”  which  raises  important  questions  about  the 
adequacy  of  the  model  relative  to  the  imaging  physics  The 
improvement  in  the  elastic  property  image  with  the  addition 
of  attenuation  effects  demonstrated  in  Fig.  6  suggests  that 
data-model  match  is  critical.  In  this  regard  it  is  important  to 
note  that  the  MR  displacement  field  is  acquired  during  the 
first  three  harmonic  cycles  of  the  displacement  stimulation 
which  is  intermittently  applied  with  the  pulse  sequencing 
(see  Fi>’.  2)  as  the  data  associated  with  each  pixel/voxel  is 
generated.  Hence,  it  is  quite  possible  that  the  measured  dis- 
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Fig.  4.  DispUicement  data  generated  by  the  phase  contrast  MR  imaging 
metliod  for  tlie  gel  sample  .shown  in  Fig.  5.  Image  A  shows  the  x  direction 
displacement  while  image  B  shows  the  y  direction  displacement,  both  in 


.Placements  do  not  reflect  fully  developed  wavefronts  which 
iiave  reached  the  dynamic  steady  state,  creating  a  conflict 
with  the  harmonic  motion  assumption  used  in  the  inversion 
model. 

To  investigate  this  possibility,  we  simulated  the  transient 
behavior  of  the  MR  data  acquisition  using  a  time-domain 
solution  of  Eq.  (4)  where  the  stimulation  was  applied  to  the 
walls  of  the  phantom  for  three  cycles  at  100  Hz.  The  time 
history  of  each  node  in  the  finite  element  mesh  which  com- 
orises  the  measured  displacement  field  was  then  Fourier 
iransformed  and  the  100-Hz  component  (magnitude  and 
phase)  was  extracted  in  order  to  represent  the  MR  displace¬ 
ment  measurements.  This  data  was  then  supplied  to  our  fre¬ 
quency  domain  inversion  algorithm  (after  being  corrupted 
with  10%  added  noise),  which  produced  the  elastic  property 
image  shown  in  image  A  of  Fig.  7.  This  image  is  strikingly 
sinfilar  in  character  to  the  images  in  Fig.  6.  In  image  B  of 
Fi".  7  we  also  show  the  case  where  the  simulated  measure¬ 
ment  data  is  taken  from  the  steady-state  solution  (with  10% 
added  noise)  and  inverted  with  the  algorithm.  The  improve¬ 
ment  in  the  recovery  of  the  elastic  property  distribution  is 
remarkable  and  suggests  that  the  degradation  of  the  image  in 
Fin.  6  results  in  part  from  data-model  mismatch  between  the 
transient  motion  which  is  measured  and  the  time-harmonic 


Fig.  5.  Stanthircl  MR  contrast  image  of  a  heterogeneous  get  sample  where 
image  B  includes  a  superimposed  outline  of  ilie  hard  inclusion. 


model  which  is  assumed.  It  is  also  interesting  to  note  that  the 
recovery  of  the  elastic  properties  is  quantitative  for  the 
steady-state  data  image  (image  B)  whereas  the  reconstructed 
properties  are  an  order  of  magnitude  too  high  for  the  tran¬ 
sient  data  image  (image  A),  presumably  due  to  some  form  of 
artificial  hardening  which  tries  to  compensate  for  the  data- 
model  mismatch.  The  reconstructions  from  experimental 
data  in  Fig.  6  exhibit  the  same  character  (i.e.,  an  effective 
hardening  of  the  embedded  inclusion). 

V.  CONCLUSIONS 

Previous  simulations^^  indicate  that  quantitative  image  re¬ 
construction  of  elastic  property  distributions  can  be  recov¬ 
ered  in  the  presence  of  considerable  measurement  noise  (up 
to  15%)  using  an  overlapping  zone  finite  element  inversion 
technique.  Reconstructed  property  resolution  can  be  pre¬ 
served  at  the  MR  measurement  level  (i.e.,  pixel  resolution) 
with  this  approach.  Initial  experience  with  reconstructions  of 
actual  MR  data  in  heterogeneous  gel  phantoms  consisting  of 
harder  inclusions  embedded  in  a  setter  background  is  en- 
couracing  and  shows  that  the  localized  increase  in  stiffness 
can  be  found  (despite  a  lack  ol  conventional  MR  image  con¬ 
trast  or  contrast  observable  directly  in  the  displacement  field 
itself).  In  this  regard,  the  incorporation  of  displacement  at¬ 
tenuation  through  the  addition  of  a  Maxwellian  viscoelastic 
term  in  the  governing  computational  model  has  been  shown 
to  improve  the  recovery  of  the  size,  shape,  and  location  of 
the  hard  inclusion.  However,  simulations  of  the  potential 
transient  effects  in  the  MR  displacement  measurements  indi¬ 
cate  that  the  data-model  mismatch  between  the  MR  data  ac¬ 
quisition  method  and  the  assumed  steady-state  image  recon- 
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Fio.  7.  Two  inversions  of  the  simulated  data  set  where  image  A  shows  the 
reconstruction  based  on  displacement  data  generated  in  three  cycles  of  mo¬ 
tion  starting  from  rest  [kPa]  while  image  B  shows  the  same  reconstruction 
process  carried  out  on  steady-state  data  [kPa]. 
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Fig.  6.  Two  inversions  of  the  MR  displacement  data  shown  in  Fig.  4,  Image 
A  shows  the  reconstruction  based  on  the  undamped  linear  elastic  model  in 
Eq.  (3)  [kPa]  and  images  B  and  C  s!u)\v  the  resulting  elasticity  (image  B) 
atid  Mawvellian  damping  (image  C)  distributions  of  the  viscoelastic  model 
[Eq.  (4)]  inversion,  in  units  of  [kPa]  and  [kg/s  lc6],  respectively. 


stniction  model  may  he  the  most  important  factor  which 
presently  limits  the  recovered  property  image  quality. 

Two  immediate  pathways  for  resolving  this  inconsistency 
present  themselves.  As  demonstrated  in  the  transient  simula¬ 
tion  reported  here,  the  MR  displacement  encoding  gradient 
imaging  process  can  be  modeled  by  taking  the  frequency 
component  of  the  transient  displacement  solution  related  to 
the  harmonic  excitation.  By  incorporating  this  process  into 


the  subzone  inversion  technique,  the  transient  effects  present 
in  displacement  images  generated  in  three  or  four  cycles  of 
actuation  could  be  modeled,  presumably  leading  to  hi^h 
quality  inversions  from  transient  data  sets  of  this  type.  Alter¬ 
natively,  by  separating  the  driving  signal  for  the  mechanical 
actuator  from  the  gradient  signal,  but  maintaining  phase  co¬ 
herency  with  the  MR  gradients  through  the  use  of  a  phase- 
lock  loop,  mechanical  excitation  could  be  applied  throughout 
the  image  acquisition  process.  This  would  assure  that  the 
resulting  displacement  images  represent  the  steady-state  mo¬ 
tion  of  the  sample,  which  should  also  lead  to  high  qua.  .> 
property  reconstructions.  Which  approach  provides  the  best 
image  reconstructions  while  still  offering  a  convenient  and 
feasible  imaging  process  remains  to  be  determined.  At  this 
point  both  of  the  strategies  mentioned  here  remain  as  viable 
options. 
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ABSTRACT:  Noise  reduction  using  monotonic  fits  between  extrema 
has  been  shown  to  work  weil  on  images,  especially  those  with  very 
iow  signai-to-noise  ratios  (SNRs).  in  this  articie  we  wili  expiore  three 
appiications  of  monotonic  noise  reduction  in  magnetic  resonance 
imaging  (MRI).  The  first  application  is  reducing  noise  In  function  MRI 
(fMRI)  studies.  Reduced  noise  ailows  greater  flexibiiity.  For  exampie, 
it  aliows  the  activated  regions  to  be  identified  using  noisier  images 
acquired  in  iess  time  or  fewer  cycles  of  stimulation.  Activation  maps 
were  calcuiated  from  the  images  after  noise  reduction  had  been 
appiied  to  each  image  in  the  series.  The  parameters  used  in  the  noise 
reduction  were  optimized  so  the  images  produced  best  matched  the 
average  of  the  entire  series.  The  CNR  was  improved  significantiy  in 
the  activation  maps.  The  resuits  can  be  extended  to  any  other  fMRi 
paradigm.  The  second  appiication  was  reducing  noise  in  complex 
data  to  improve  the  SNR  of  the  phase  in  the  compiex  MRI  image.  The 
error  in  the  phase  was  reduced  by  a  factor  of  three  in  the  simuiations 
shown,  in  the  third  appiication,  we  introduced  a  simple  contrast- 
enhancement  method  using  monotonic  noise  reduction.  To  enhance 
contrast,  the  coarse  features  were  reduced  in  size;  the  smaiier  size 
features  were  increased  in  size;  very  smaii  features  that  are  likeiy  to  be 
noise  were  attenuated.  The  result  is  a  simpie,  effective  method  of 
improving  the  contrast  of  features  of  a  seiected  size  in  images  with  no 
false  features  introduced.  ©  1 999  John  Wiley  &  Sons,  Inc.  Int  J  Imaging  Syst 
Technol,  10,  177-185,  1999 

Key  words:  noise  reduction;  fMRi;  contrast  enhancement;  phase 
estimation 


I.  INTRODUCTION 

Monotonic  noise  reduction  is  a  relatively  new  method  that  is  effec¬ 
tive  at  lower  signal-to-noise  ratios  (SNRs)  than  other  methods 
(Weaver,  1997).  It  also  avoids  Gibb’s  ringing  and  blurring,  two 
common  problems  that  afflict  currently  used  methods.  Noise  reduc¬ 
tion  filters  and  denoising  methods  fit  many  adjacent  data  points  to  a 
set  of  basis  functions.  This  change  of  basis  allows  the  parts  of  the 
data  most  likely  to  be  noise  to  be  suppressed.  The  effectiveness  of 
the  particular  method  depends  on  how  well  the  basis  separates  noise 
from  data.  In  high  noise  regions,  projection  onto  the  basis  functions 
that  average  large  numbers  of  pixels  represents  the  data  most  accu¬ 
rately  and  therefore  has  the  least  noise.  Representing  the  data  with 
only  those  basis  functions  leads  directly  to  bandlimiting  loss  of 
spatial  resolution  and  undersampling  artifacts.  Monotonic  noise 
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reduction  does  not  fit  the  data  with  a  limited  set  of  basis  functions, 
so  no  undersampling  or  ringing  artifacts  are  introduced  and  there  is 
no  bandlimiting.  It  simply  constrains  the  data  to  be  monotonic 
between  extrema. 

A.  Monotonic  Noise  Reduction  in  One  Dimension.  Mono¬ 
tonic  noise  reduction  is  simplest  in  one  dimension.  An  elegant 
algorithm  (Demetriou,  1990;  Demetriou  and  Powell,  1991)  that 
produces  the  monotonically  increasing  series  that  fits  a  data  set  best 
in  a  least-squares  sense  is  used.  The  algorithm  travels  through  the 
data  until  it  finds  a  decrease.  A  decrease  violates  the  monotonic 
condition.  When  it  finds  a  decrease,  the  algorithm  averages  that 
decreasing  data  point  with  as  many  previous  points  as  is  necessary 
to  establish  monotonic  increase.  The  algorithm  then  continues  the 
process  by  searching  for  the  next  decrease.  When  all  of  the  decreases 
are  averaged  away,  the  algorithm  is  completed. 

The  primary  advantage  of  this  algorithm  for  noise  reduction  is 
that  it  smoothes  the  data  as  little  as  possible  and  does  not  blur  edges. 
The  algorithm  leaves  all  increases  unchanged;  both  sharp  and 
smooth  increases  remain  unchanged  so  no  smoothing  occurs  at  all. 
The  algorithm  averages  just  as  many  data  points  as  necessary  to  get 
a  monotonically  increasing  series  so  the  data  are  smoothed  as  little 
as  possible. 

However,  not  all  signals  are  monotonically  increasing  functions. 
A  monotonically  decreasing  series  can  be  easily  obtained  with  the 
same  algorithm  by  reversing  the  order  of  the  data  or  by  subtracting 
the  data  from  some  large  constant.  The  algorithm  can  also  be  used 
to  produce  a  piecewise  monotonic  series  by  breaking  the  data  into 
regions  between  extrema  that  are  each  fit  to  increasing  or  decreasing 
monotonic  series. 

The  key  to  fitting  piecewise  monotonic  series  is  to  identify  the 
extrema  correctly.  In  this  sense,  monotonic  noise  reduction  is  very 
similar  to  wavelet  denoising.  The  key  in  both  cases  is  in  identifying 
the  extrema.  In  wavelet  denoising,  extrema  in  the  gradient  are 
required  and  in  monotonic  noise  reduction  the  signal  peaks  are 
required.  To  select  significant  extrema,  we  have  used  a  Fourier  filter 
to  smooth  small  extrema  away.  Then  only  extrema  that  are  different 
from  previous  and  subsequent  extrema  by  more  than  a  user-defined 
threshold  are  used.  There  are  two  parameters  that  the  user  must  set 
to  find  the  extrema:  the  bandwidth  allowed  through  the  Fourier  filter 
and  the  minimum  difference  between  extrema  that  is  allowed.  When 
the  extrema  are  found,  noise  can  be  suppressed  by  forcing  mono¬ 
tonic  change  between  the  extrema.  The  algorithm  in  one  dimension 
is  shown  in  Figure  1. 
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Figure  1.  The  one-dimensional  monotonic  noise-reduction  algo¬ 
rithm.  (Top  to  bottom)  The  original  signal,  the  signal  plus  added  noise, 
the  blurred  signal  with  the  selected  extrema  shown  as  vertical  lines, 
and  the  results  of  monotonic  noise  reduction.  The  SNR  is  10;  the 
standard  deviation  of  the  noise  is  one  tenth  and  the  peak  signal  is  one 
so  the  total  signal  energy  is  less  than  1 0  times  the  total  noise  energy. 


B.  Monotonic  Noise  Reduction  in  Two  Dimensions,  fhe 

extension  of  the  algorithm  to  multiple  dimensions  is  accomplished 
by  ordering  the  pixels  in  the  image  into  a  one-dimcnsional  series 
where  the  one-dimcnsional  algorithm  can  be  used.  In  previous 
versions  of  this  algorithm,  the  one-dimensional  algorithm  was  run 
on  the  rows  and  then  the  columns  of  rotated  versions  of  the  image 
(Weaver,  1997b).  As  many  rotations  as  needed  were  used.  However, 
most  of  the  noise  is  eliminated  using  just  four  angles.  Therefore,  a 
faster  version  of  the  algorithm  that  uses  the  one-dimcnsional  algo¬ 
rithm  on  the  rows,  columns,  and  diagonals  of  the  image  is  used  here. 
Figures  2  and  3  illustrate  the  steps  used  in  the  algorithm.  No 
interpolation  is  required  to  rotate  the  image  to  other  angles,  which 
speeds  up  the  algorithm  significantly.  This  version  of  the  algorithm 
was  used  in  all  of  the  applications  presented  here.  This  algorithm 
was  coded  in  Matlab@  with  two  subroutines  coded  in  C.  It  runs  in 
10-20  s  on  a  256  X  256-pixel  image  on  a  DEC  Alpha @.  Earlier 
versions  of  the  algorithm  required  2-3  min.  In  both  cases,  the  time 
required  depends  on  the  number  of  extrema  identified.  More  extrema 
increases  the  run  time.  Interactive  languages  like  Matlab@  are 
significantly  slower  than  Fortran  or  C,  so  execution  times  could  be 
reduced  significantly. 


Figure  2.  The  selection  of  extrema  in  the  two-dimensional  monotonic  noise  reduction  algorithm.  The  top  image  is  the  original  image  with 
significant  noise.  The  Fourier  blurred  image  is  shown  below  the  original  image.  The  extrema  are  selected  along  the  columns,  rows,  and  both 
diagonals  of  the  two-dimensionally  blurred  image.  The  four  images  at  the  bottom  show  the  position  of  the  selected  extrema  for  all  four  directions. 
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Figure  3.  Application  of  the  one-dimensional  algorithm  in  the  two 
dimensional  monotonic  noise-reduction  algorithm.  The  pixels  be¬ 
tween  extrema  along  the  columns  are  made  monotonic.  The  rows  of 
the  processed  image  are  then  processed  using  the  one-dimensional 
algorithm.  This  forms  one  intermediate  image.  Processing  the  rows 
and  then  the  columns  forms  another  intermediate  image.  Processing 
the  upper  left  to  lower  right  diagonal  followed  by  the  lower  left  to 
upper  right  diagonal  forms  another  intermediate  images.  Processing 
the  lower  left  to  upper  right  diagonal  followed  by  the  upper  left  to 
lower  right  diagonal  forms  the  last  intermediate  image.  All  four  inter¬ 
mediate  images  are  averaged  together  to  form  the  final  processed 
image. 


II.  REDUCING  NOISE  IN  FMRI  ACTIVATION  MAPS 

In  functional  magnetic  resonance  imaging  (fMRI),  multiple  serial 
images  of  the  same  anatomy  are  obtained  during  a  series  of  mental 
tasks  (Frank  et  al.,  1998;  Jezzard  and  Song,  1996;  Turner,  1997; 
Weisskoff,  1995).  Small  changes  in  the  signal  of  brain  tissue  reflect 
the  small  blood  flow  changes  and  oxygenation  levels  that  are  asso¬ 
ciated  with  activation.  At  1.5  T  the  changes  in  signal  are  around  3% 
(Turner  et  al.,  1993),  which  is  the  same  size  as  the  random  noise  in 
good-quality  images.  Therefore,  only  changes  that  are  correlated 
with  the  mental  tasks  over  many  task  cycles  signify  activation.  Noise 
in  the  images  limits  fMRI  studies  in  several  ways.  The  images  are 
generally  low  resolution  to  obtain  sufficient  SNR  to  detect  the 
activations.  As  much  imaging  time  as  possible  is  used  to  obtain  the 
required  SNR.  Reducing  the  imaging  time  would  reduce  the  SNR, 
but  it  would  mitigate  several  other  problems.  It  might  be  possible  to 


reduce  artifacts  and  it  would  allow  more  flexibility  in  task  design. 
Many  mental  tasks  are  difficult  to  turn  on  and  off  consistently  over 
several  cycles.  For  example,  memory  tasks  are  difficult  to  turn  off 
because  it  is  difficult  to  stop  thinking  about  something.  Distraction 
can  be  effective  but  is  not  consistent  and  distraction  introduces  other 
activations.  On  the  other  hand,  many  activations  are  difficult  to 
maintain  in  a  controlled  manner.  Concentration  can  wander  to  the 
noise  of  the  scanner  or  other  unintended  subjects  very  easily.  There¬ 
fore,  reducing  the  sensitivity  to  noise  is  important.  We  have  used 
monotonic  noise  reduction  to  eliminate  some  of  the  random  noise  in 
the  images  so  the  activation-related  signal  changes  are  more  easily 
identified  from  noisier  images  that  can  take  less  time  to  acquire. 

A.  Methods.  We  used  the  monotonic  noise  reduction  algorithm 
described  above  on  each  image  in  an  fMRI  motor  activation  study. 
The  unusual  aspect  of  filtering  fMRI  images  is  that  they  are  all 
nearly  identical  so  the  average  image  is  an  excellent,  low-noise 
estimate  of  each  individual  image.  We  used  the  average  image  to 
estimate  the  best  parameters  to  use  in  the  noise  reduction  algorithm. 
The  parameters  that  made  one  of  the  images  in  the  middle  of  the 
series  most  like  the  average  image  were  used.  The  Nelder-Mead- 
type  simplex  search  method  implemented  in  Matlab@  was  used  to 
find  the  optimum  parameters  for  the  monotonic  noise-reduction 
method.  Then,  each  image  in  the  series  was  processed  using  the 
same  parameters.  Simple  z -scores  of  the  images  were  used  to  form 
the  activation  map. 

We  used  a  simple  motor  activation  study  that  clearly  showed 
both  right-  and  left-hand  activation.  We  then  added  normally  dis¬ 
tributed  random  noise  to  the  images.  The  activation  maps  were 
calculated  with  and  without  monotonic  noise  reduction.  The  average 
z-scores  in  the  activated  regions  and  in  the  background  brain  were 
calculated  to  evaluate  how  well  the  activation  could  be  identified. 
The  activated  regions  were  those  that  had  a  z-score  above  3  in  the 
original  activation  maps.  The  background  region  was  the  entire  brain 
except  rectangular  regions  centered  on  the  right  and  left  activations. 
The  activation  maps  were  calculated  with  nine  amounts  of  added 
noise.  The  standard  deviation  of  the  noise  ranged  from  0%  to  25% 
of  the  maximum  signal  in  the  image.  The  noise  energy  added  ranged 
from  0%  of  the  total  image  energy  to  93%  of  the  total  image  energy. 

The  motor  study  we  used  is  a  simple  one  that  is  used  in  a 
laboratory  exercise  for  engineering  students.  Although  the  study 
design  is  old,  the  random  noise  is  representative  of  almost  any  fMRI 
study.  Therefore,  the  noise  reduction  shown  using  this  study  is 
representative  of  noise  reduction  using  other  fMRI  studies.  The 
study  was  performed  on  a  1.5-T  General  Electric  Signa  system  with 
10  mT/m  gradients.  A  single  axial  plane  4  cm  below  the  apex  of  the 
cranium  was  imaged.  Six  cycles  of  rest  and  motor  activation  were 
obtained.  Each  cycle  of  images  consisted  of  five  resting  images 
followed  by  five  images  taken  during  tapping  of  the  right  hand 
followed  by  five  images  taken  during  tapping  of  the  left  hand.  Each 
image  was  a  256  X  256,  24-cm  field  of  view,  gradient-echo  image 
with  a  TR  of  70  ms  and  TE  of  40  ms.  The  images  were  strongly  T2* 
weighted,  z -scores  for  each  pixel  were  used  to  identify  activated 
regions.  The  standard  deviation  of  the  90  images  was  used  as  a 
measure  of  the  total  variation  in  signal.  To  obtain  the  signal  corre¬ 
lated  with  right  hand  finger  tapping,  the  average  of  the  rest  and  left 
hand  tapping  was  subtracted  from  the  average  of  the  right  hand 
finger  tapping.  The  activation  map  for  the  right  hand  motor  task  was 
the  ratio  of  the  signal  correlated  with  right  hand  finger  tapping  over 
the  total  variation  in  signal.  Similarly,  the  activation  map  for  the  left 
hand  motor  task  was  the  ratio  of  the  signal  correlated  with  left  hand 
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Figure  4.  The  activation  maps  for  left  and  right  motor  activity  obtained  from  six  cycles  of  a  simple  motor  task.  The  maps  in  each  column  were 
obtained  from  the  original  data  with  different  amounts  of  added  noise.  Ten  percent  added  noise  means  that  the  standard  deviation  of  the  random 
noise  is  10%  of  the  maximum  pixel  intensity  in  the  average  image.  The  top  row  are  the  activation  maps  for  the  right  hand  task.  The  second  row 
are  the  activation  maps  for  the  left  hand  task.  The  third  and  fourth  rows  are  the  activation  maps  for  the  right  hand  and  left  hand,  respectively, 
when  monotonic  noise  reduction  is  applied  to  the  raw  images  before  the  z-scores  are  calculated.  False  activations  appear  in  both  right  and  left 
maps  with  10%  noise  when  noise  reduction  is  not  used.  Monotonic  noise  reduction  allows  the  activations  to  be  identified  when  as  much  as  20% 
added  noise;  in  this  case,  the  noise  energy  is  essentially  half  of  the  total  energy  in  the  image.  The  regions  identified  were  those  with  at  least  one 
pixel  with  the  maximum  value  in  the  true  activated  region.  All  pixels  in  the  region  were  above  the  mean  z-score  in  the  true  activated  region.  This 
method  ensured  that  the  true  activated  region  was  identified,  but  it  shows  how  unique  that  region  is. 


finger  tapping  over  the  total  variation  in  signal.  The  images  were 
aligned  using  a  multiscale,  rigid  alignment  method  developed  by 
Unser  et  al.  (1993)  and  Kostelec  et  al.  (1998).  Alignment  did  not 
change  the  results  dramatically. 

B.  Results.  The  z-score  maps  calculated  from  the  original  images 
were  used  as  the  true  position  of  the  motor  strips.  The  first  column 
of  Figure  4  shows  the  activation  maps  obtained  from  all  six  cycles 
of  the  right-  and  left-hand  motor  task  both  with  and  without  noise 
reduction.  When  the  standard  deviation  of  the  added  noise  was  10% 
of  the  maximum  signal  value,  spurious  activations  were  seen  in  the 
right  and  left  unfiltered  maps.  There  were  no  spurious  activations 
when  noise  reduction  was  used.  When  the  standard  deviation  of  the 
added  noise  was  20%  of  the  maximum  signal  value,  the  spurious 
activations  eliminated  all  possibility  of  identifying  the  correct  acti¬ 


vated  regions  when  noise  reduction  was  not  used.  There  were  no 
spurious  activations  in  the  right-hand  activation  map  when  noise 
reduction  was  used,  and  few  in  the  left.  The  size  of  the  activated 
regions  increased  when  noise  reduction  was  used  because  the  cutoffs 
used  to  identify  the  activation  were  set  by  the  unprocessed  maps. 

Figures  5  and  6  show  the  mean  z-scores  in  the  activated  regions 
and  in  the  background  for  right-  and  left-hand  tasks,  respectively, 
using  monotonic  noise  reduction  increased  the  z-score  in  the  acti¬ 
vated  region  by  as  much  as  128%  for  the  right-hand  task  and  143%^ 
for  the  left-hand  task.  The  average  increase  was  66%  for  the  right- 
hand  task  and  75%  for  the  left-hand  task.  The  average  z-scores  in  the 
background  region  remained  essentially  unchanged. 

Noise  reduction  algorithms  were  only  effective  for  random  noise. 
Systematic  sources  of  noise  such  as  patient  motion  or  the  cerebro¬ 
spinal  fluid  (CSF)  moving  cannot  be  reduced  with  these  methods. 
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Figure  5.  The  mean  z-scores  of  pixels  in  the  region  activated  during 
right-hand  activity  as  a  function  of  added  noise.  (Top  curve)  When 
noise  reduction  is  used.  (Middle  curve)  When  the  original  images  were 
used.  The  two  curves  along  the  bottom  are  the  mean  z-scores  in  the 
background  region.  The  background  did  not  change  significantly.  The 
background  Is  the  entire  brain  excluding  a  rectangular  region  around 
the  right  and  left  hand  activated  regions.  Noise  reduction  increased 
the  mean  z-score  by  as  much  as  128%. 


However,  the  methods  of  image  acquisition  might  be  changed  to 
speed  up  imaging  or  synchronize  it  with  the  other  sources  of  phys¬ 
iologic  noise  even  at  the  expense  of  increased  random  noise  if  noise 
reduction  is  employed.  Random  noise  reduction  is  a  useful  tool  that 
can  be  used  to  improve  fMRI  studies  when  the  random  noise  masks 
the  real  activations. 

III.  PHASE  RECOVERY 

There  are  many  applications  where  the  phase  of  the  MR  signal 
encodes  important  information.  For  example,  the  phase  can  be  used 
to  generate  frequency  maps.  Frequency  maps  might  be  useful  in 
fMRI  to  separate  the  effects  of  large-scale  susceptibility  changes 
caused  by  larger  vessels  from  intravoxel  susceptibility  changes 
caused  by  changes  in  capillary  blood  flow.  The  phase  can  also  be 
used  to  encode  velocity  information  in  phase-contrast  imaging  (Mo¬ 
ran,  1982;  Van  Dijk,  1984).  In  elastography  (Muthupillai  and  Eh- 
man,  1996;  Plewes  et  al.,  1995;  Chenevert  et  al.,  1998),  the  phase 
encodes  the  displacement  of  tissue  produced  by  low-frequency 
vibrations.  In  many  of  these  cases,  these  measurements  have  rela¬ 
tively  low  SNRs  because  they  must  be  obtained  with  fast  imaging 
methods. 

Both  the  real  and  imaginary  parts  of  the  complex  image  are 
distorted  by  additive  normally  distributed  random  noise.  The  noise 
becomes  Rician  only  after  the  magnitude  is  taken  (Macovski,  1996). 
The  real  and  imaginary  parts  of  the  signal  can  be  filtered  separately 
with  no  added  bias  because  the  noise  is  independent.  The  uncertainty 
in  the  phase  angle  at  each  pixel  still  depends  on  the  magnitude  of  the 
signal,  but  reduced  noise  decreases  that  uncertainty. 

It  is  not  uncommon  to  repetitively  scan  of  the  same  physical 
position.  In  elasticity  measurements  using  phase-contrast  estimates 
of  periodic  motion,  three  scans  can  be  taken  to  measure  displace¬ 
ment  in  all  three  directions.  The  magnitude  of  the  image  is  identical 
in  all  three  scans;  only  the  phase  changes.  Again  as  in  the  fMRI 


experiment,  the  average  provides  a  good  estimate  of  each  image  so 
we  can  use  the  same  method  of  finding  the  optimum  parameters  for 
the  monotonic  noise  reduction  algorithm. 

A.  Methods.  Noise  in  the  real  and  imaginary  parts  of  the  signal  is 
suppressed  separately  using  monotonic  noise  reduction.  The  phase  is 
calculated  from  the  resulting  complex  signal. 

A  one-dimensional  simulation  has  been  used  to  demonstrate  the 
potential.  The  phase  varied  linearly  from  tt  to  -tt  across  the  signal. 
The  total  signal  energy  over  the  total  noise  energy  is  0.38.  The 
standard  deviation  of  the  noise  in  both  the  real  and  imaginary 
channels  is  0.25  and  the  peak  signal  is  one. 

For  the  two-dimensional  case,  it  was  assumed  that  the  magnitude 
of  the  image  was  known.  In  that  case,  the  average  can  be  used  to  find 
the  optimum  parameters  in  the  filtration.  The  threshold  and  band¬ 
width  used  to  find  the  extrema  were  selected  to  minimize  the 
difference  between  the  magnitude  of  the  filtered  image  and  the 
magnitude  of  the  average  image.  The  Nelder-Mead  type  simplex 
search  method  implemented  in  Matlab@  was  used  to  find  the 
optimum  parameters  for  the  monotonic  noise  reduction  method.  For 
demonstration,  we  assumed  that  the  magnitude  was  a  clinical  axial 
image.  The  phase  was  simulated  with  a  quadratic  polynomial.  Nor¬ 
mally  distributed  random  noise  was  added  to  both  images  to  disturb 
the  phase  and  magnitude  of  the  image. 

B.  Results.  The  one-dimensional  signal  is  shown  in  Figure  7.  The 
standard  deviation  of  the  error  in  phase  was  reduced  from  3 1  degrees 
in  the  original  noisy  data  to  seven  degrees  following  monotonic 
noise  reduction. 

The  two  dimensional  results  also  depend  on  the  initial  SNR  in  the 
complex  data.  Figure  8  shows  the  phase  of  the  original  image,  the 
noisy  image  and  the  results  of  noise  reduction.  The  mean  error  in  the 
phase  was  decreased  by  a  factor  of  almost  three  from  58  degrees  to 
21  degrees. 

IV.  CONTRAST  ENHANCEMENT 

Contrast  enhancement  can  be  extremely  useful  in  a  wide  variety  of 
applications.  It  can  be  used  to  improve  the  conspicuity  of  low- 


Figure  6.  The  mean  z-scores  of  pixels  in  the  region  activated  during 
left-hand  activity  as  a  function  of  added  noise.  The  order  of  curves  is 
the  same  as  in  Figure  5.  Noise  reduction  increased  the  mean  z-score 
by  as  much  as  143%. 
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Figure  7.  (a)  The  real  and  imaginary  parts  of  the  simulated  signal,  (b)  The  signal  with  added  noise,  (c)  The  signal  recovered  from  the  noisy  data 
using  monotonic  noise  reduction,  (d)  The  top  plot  is  the  phase  of  the  unprocessed  signal  plotted  with  the  true  phase.  The  bottom  shows  the 
phase  recovered  from  monotonic  noise  reduction  compared  to  the  original  phase  is  shown  on  the  bottom.  The  noise  is  very  accurately  recovered 
in  regions  where  the  signal  is  large  enough.  In  regions  where  the  signal  is  almost  zero,  the  phase  is  not  accurate. 


contrast  features  making  them  harder  to  overlook.  Contrast  enhance¬ 
ment  has  been  used  in  mammography  to  make  lesions  more  prom¬ 
inent  (Lu  et  al„  1994;  Li  et  al.,  1997).  Contrast  enhancement  can  be 
used  on  digitized  chest  images  (Rehm  et  al.,  1990;  Souto  et  al., 
1992;  McNitt-Gray  et  al.,  1993)  and  routine  bone  images  (Ogoda  et 
al.,  1997),  or  it  could  be  used  in  radiation  therapy  to  improve  the 
contrast  of  port  films. 

Contrast  enhancement  can  be  accomplished  with  histogram  mod¬ 
ification  (Gonzalez  and  Wintz,  1977)  or  edge-based  techniques 
(Beghdadi  and  LeNegrate,  1989;  Neycenssac,  1993).  Histogram 
equalization  performs  poorly  with  large  uniform  regions.  Several 
edge-based  methods  (Barrett  and  Swindell,  1981)  can  be  reduced  to 
variations  of  the  blurred-mask  algorithm  (Bednarek  and  Rubin. 
1991;  Jain,  1989).  Blurred-mask  or  unsharp-mask  methods  can  be 
used  to  amplify  different  size  features  by  different  amounts.  Other 


techniques  that  use  the  wavelet  transform  can  amplify  features  of 
different  sizes  by  different  amounts  as  well  (Lu  et  al.,  1994b).  We 
have  used  monotonic  noise  suppression  in  a  blurred-mask  method  to 
enhance  contrast.  The  primary  advantages  of  the  monotonic  method 
are  that  it  works  well  at  lower  SNRs,  is  not  bandlimited,  and  does 
not  introduce  ringing  in  the  image.  Ringing  produced  by  undcrsam- 
pling  can  be  amplified  and  introduce  false  lesions.  Unsharp-mask 
methods  also  allow  the  basic  intensity  distribution  to  remain  similar, 
so  the  look  and  feel  of  the  contrast-enhanced  image  are  the  same  as 
those  of  the  original  image. 

A.  Methods.  We  used  monotonic  noise  reduction  in  a  variation  of 
the  blurred-mask  method  to  enhance  contrast.  Two  filtered  images, 
I,  and  L,  are  obtained  using  the  standard  noise-reduction  algorithm. 
The  more  heavily  filtered  image,  L,  is  used  to  reduce  the  baseline 
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Figure  8.  The  phase  of  the  noiseless  complex  Image  is  the  upper 
right  image.  The  other  two  images  in  the  top  row  are  the  phase  of  the 
noisy  image  and  phase  error  in  the  noisy  image.  The  two  images  in  the 
bottom  row  are  the  phase  of  the  monotonically  noise  suppressed 
image  and  the  phase  error.  The  noise  in  the  monotonically  noise 
suppressed  image  is  reduced  significantly.  The  error  in  the  phase  is 
improved  by  a  factor  of  almost  three,  from  58  degrees  In  the  noisy 
image  to  21  degrees  in  the  monotonically  noise  suppressed  image. 


variations  in  the  image.  The  details  to  be  enhanced  are  left  in  the  first 
image,  Ij,  and  filtered  out  of  the  second  image,  I2.  Those  details  are 
found  by  subtracting  I2  from  I^.  However,  the  general  brightness 
pattern  of  the  image  needs  to  be  maintained  so  the  details  are 
amplified  and  added  to 

Ice  =  II  +  “  I2) 


Figure  9.  The  top  curve  is  the  phase  of  the  128th  row  of  the 
noiseless  image  shown  in  Figure  8.  The  middle  row  is  the  phase  of  the 
128th  row  of  the  noisy  image.  The  bottom  row  is  the  phase  of  the 
128th  row  of  the  monotonically  noise  suppressed  Image. 


of  filtered  images  obtained  using  increasingly  large  thresholds  so 
each  image  has  fewer  small  features  than  previous  ones.  Subtracting 
pairs  of  adjacent  images  can  separate  the  features  of  each  size  range. 
Each  subtracted  pair  can  be  amplified  by  any  desired  amount  in  the 
contrast-enhanced  image: 


where  a  is  the  factor  used  to  vary  the  contrast  enhancement.  Factors 

of  10-15  are  commonly  used  in  our  experience.  «-i 

Features  of  different  sizes  can  be  amplified  by  different  amounts  =  Ii  F  ^ 

if  multiple  filtered  images  are  obtained.  Let  I2, . . . ,  I„  be  a  series 


Figure  10,  Original  chest  image  is  shown  on  the  left  and  the  contrast  enhanced  image  is  shown  on  the  right.  The  image  was  enhanced  to 
visualize  the  vessels  in  the  lung  and  the  pacemaker  leads.  The  vessels  in  the  lung  and  the  sharp  outlines  of  the  vertebrae  In  the  mediastinum 
are  all  much  more  clearly  visible  in  the  contrast-enhanced  image.  The  leads  are  visible  in  the  lung  and  as  they  cross  the  mediastinum. 
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B.  Results.  We  enhanced  the  contrast  in  a  digitized  chest  image 
to  improve  the  conspicuity  of  the  vessels  in  the  lung  and  the 
pacemaker  lines.  The  original  and  contrast  enhanced  images  are 
shown  in  Figure  10.  It  is  clear  that  the  vessels  are  well  depicted  and 
there  are  no  spurious  signals.  The  vertebrae  and  the  ribs  are  also  seen 
much  more  clearly  in  the  contrast-enhanced  image. 

We  also  enhanced  the  contrast  in  a  mammogram.  Both  the 
original  image  and  the  contrast  enhanced  image  are  shown  in  Figure 

11.  The  vessels  throughout  the  breast  are  much  more  clearly  de¬ 
picted.  A  plot  of  the  150th  column  of  both  images  is  shown  in  Figure 

12.  The  plot  allows  the  reader  to  match  features  in  the  original  image 
with  features  in  the  contrast-enhanced  image.  No  extra  peaks  were 
introduced  by  the  contrast  enhancement,  and  the  general  intensity 
distribution  of  the  image  was  maintained. 

V.  CONCLUSIONS 

Monotonic  noise  reduction  is  useful  in  situations  where  noise  is 
large,  and  it  is  especially  useful  when  significant  information  about 
the  extrema  exists.  We  have  shown  three  examples  of  the  utility  of 
monotonic  noise  reduction.  Suppressing  noise  prior  to  calculating 
the  z-scores  significantly  enhances  fMRI  activation  maps.  The  CNR 
can  be  more  than  doubled.  The  average  image  provides  an  effective 
way  to  estimate  the  optimum  parameters  for  the  noise-reduction 
algorithm.  The  phase  error  can  be  improved  significantly  by  filtering 
both  real  and  imaginary  parts  of  the  reconstructed  image  before 
calculating  the  phase.  If  only  the  phase  changes  between  multiple 
measurements  of  the  same  physical  region,  the  accurate  estimate  of 
the  magnitude  of  the  image  can  be  used  to  improve  the  estimate  of 
the  phase  even  more.  Finally,  a  simple  robust  method  of  contrast 
enhancement  using  monotonic  noise  reduction  was  shown  to  be 
effective  in  maintaining  the  general  appearance  of  the  image  while 
improving  contrast  significantly.  There  are  many  improvements  that 
can  be  made  to  the  algorithm.  Methods  of  reducing  the  noise  leakage 


Figure  11.  A  mammogram  is  shown  on  the  left  and  the  contrast 
enhanced  mammogram  is  shown  on  the  right.  The  vessels  and  ducts 
are  much  more  clearly  delineated. 


Figure  12.  Plot  comparing  of  the  150th  column  of  the  original  and 
contrast  enhanced  breast  images.  (A)  The  original  image.  (B)  The 
contrast  enhanced  image.  (C)  Rows  100-200  of  the  original  image  in 
greater  detail.  (D)  Rows  100-200  of  the  contrast  enhanced  image  for 
comparison.  The  general  intensity  distribution  is  maintained.  The 
expanded  views  show  that  each  peak  in  the  contrast  enhanced  image 
corresponds  to  a  peak  in  the  original  image.  There  are  no  new  peaks 
and  the  position  of  the  peaks  remains  the  same.  No  spurious  features 
are  introduced  and  there  is  no  ringing  to  interfere  with  the  vessels. 

at  boundaries  are  particularly  promising.  There  are  also  many  vari¬ 
ations  of  the  basic  method  that  can  be  made  to  adapt  it  to  specific 
situations. 
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Abstract 


A  novel  method  of  reducing  noise  in  images  is  presented.  The  significant  extrema  (maxima  and 
minima)  in  the  image  are  selected  using  a  simple  low  pass  Fourier  filter.  The  method  forces  the  pixel 
values  in  the  image  to  vary  monotonically  between  the  selected  extrema.  For  example,  the  pixel  values 
in  the  filtered  image  should  decrease  monotonically  in  all  directions  from  an  isolated  maximum.  Because 
the  algorithm  that  performs  the  monotonic  fits  is  one  dimensional,  we  approximate  monotonic  change  in 
all  directions  by  doing  monotonic  fits  along  line  segments  throughout  the  image.  The  filtering  operation 
on  each  line  segment  replaces  the  pixel  values  on  that  segment  with  a  monotonic  sequence  that  fits  the 
original  pixel  values  best  in  a  least  squares  sense.  Monotonic  change  is  enforced  along  line  segments  in  as 
many  directions  as  desired.  The  method  is  simple,  reasonably  fast  and  quite  stable.  Good  results  can  be 
obtained  for  images  with  SNR’s  as  low  as  0.5. 
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L  Background 

There  are  many  methods  of  removing  noise  from  images.  The  basic  idea  is  almost  always  to  represent 
the  signal  and  noise  in  a  basis  that  separates  the  true  signal  from  the  noise  as  completely  as  possible.  In 
Fourier  filters,  the  noise  is  assumed  to  be  dominant  in  the  high  frequencies  and  the  signal  dominant  in 
the  low  frequencies  so  removing  high  frequencies  removes  mostly  noise.  However,  removing  or  damping 
the  high  frequencies  produces  bandlimited  blurring.  Wavelet  denoising  [1,  2,  3,  4,  5,  6]  adds  another 
twist  by  allowing  the  bandwidth  passed  to  be  different  at  different  positions  in  the  image.  The  high 
frequencies  are  allowed  to  pass  where  edges  have  been  identified  or  are  likely  and  the  high  frequencies  are 
more  strongly  suppressed  where  edges  are  not  likely.  Edges  produce  extrema,  in  the  wavelet  transform 
domain  so  a  wider  bandwidth  is  passed  around  large  values  of  the  wavelet  transform.  One  of  the  keys 
in  wavelet  denoising  is  deciding  when  an  edge  is  present  and  when  the  large  wavelet  response  is  from 
noise.  An  alternative  method  of  preferentially  suppressing  noise  is  presented  here.  Instead  of  identifying 
extrema  of  the  derivative  as  in  wavelet  denoising,  the  extrema  of  the  signal  itself  (peaks)  are  found;  the 
signal  is  then  forced  to  change  monotonically  between  the  extrema.  The  method  should  be  more  robust 
than  wavelet  denoising  because  averages  are  more  stable  than  differences.  However,  identifying  the  correct 
extrema  is  again  the  key  to  the  performance  of  the  method.  The  filter  is  made  possible  by  an  elegant 
algorithm  due  to  Demetriou  and  Powell  that  finds  the  monotonic  series  that  fits  data  points  best  in  a  least 
squares  sense  [7].  It  basically  averages  large  enough  groups  of  adjacent  points  in  the  series  to  achieve  a 
monotonic  progression.  Monotonically  decreasing  series  can  be  found  by  reversing  the  order  of  the  data. 
It  is  a  robust  and  relatively  fast  method.  Our  noise  reduction  method  for  images  has  two  parts:  selecting 
the  important  extrema  and  forcing  monotonic  change  between  those  extrema.  Extrema  selection  is  best 
done  in  two  dimensions  because  extra  stability  is  gained  by  smoothing  in  both  dimensions  rather  than 
only  one.  The  process  of  forcing  monotonic  change  in  all  directions  is  accomplished  by  forcing  monotonic 
change  along  many  line  segments.  We  will  first  look  at  the  relevant  properties  of  the  algorithm  that  forces 
monotonic  change  along  a  single  line  segment.  Then  we  will  describe  the  two  dimensional  algorithm  in 
detail,  show  some  results,  give  some  applications,  and  propose  future  improvements. 

II.  One  Dimensional  Examples 

The  values  of  a  function  at  points  between  two  adjacent  extrema  must  vary  monotonically  essentially 
by  definition.  Therefore,  a  function  of  one  variable  in  additive  noise  can  be  retrieved  by  identifying  the 
positions  of  the  extrema  and  fitting  the  points  bounded  by  each  pair  of  adjacent  extrema  to  a  mono  tonic 
function.  For  any  segment  bounded  by  a  pair  of  adjacent  extrema,  the  monotonic  series  that  fits  the 
points  best  in  a  least  squares  sense  can  be  found  using  the  Demetriou  and  Powell  algorithm  [7].  This 
algorithm  has  several  very  interesting  properties.  The  result  is  not  just  a  fit  to  a  limited  set  of  basis 
functions.  Any  monotonic  series  can  be  recovered.  The  result  is  not  bandlimited;  edges  are  not  blurred. 
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The  algorithm  takes  the  first  element  of  the  series  as  the  starting  value  and  advances  through  the  series 
element  by  element.  If  an  element  in  the  series  is  smaller  than  the  previous  element,  the  algorithm  averages 
past  points  in  the  series  to  bring  the  past  values  down  and  the  present  element  up  enough  to  maintain 
monotonic  increase.  The  algorithm  is  fast,  works  in  place  so  memory  requirements  are  minimal  and  is 
relatively  robust.  Figures  1  to  4  demonstrates  the  important  properties.  Figure  1  shows  a  simulated  signal 
made  up  of  three  boxcars  with  moderate  added  noise.  The  local  SNR  over  the  boxcars  is  approximately 
ten.  The  ratio  of  the  total  signal  energy  to  the  total  noise  energy  is  six.  This  is  generally  a  relatively 
poor  SNR  for  medical  imaging;  a  good,  crisp  MR  image  generally  has  a  SNR  of  at  least  80.  The  blurred 
signal  and  the  extrema  selected  are  shown  as  well  as  the  final  noise  suppressed  signal.  Two  features  of  the 
noise  suppressed  signal  should  be  noted;  the  edges  are  not  blurred  at  all  and  there  are  spurious  peaks  at 
the  boundaries  between  monotonic  segments.  Demetriou  and  Powell’s  algorithm  is  very  happy  with  sharp 
edges;  they  are  not  blurred  at  all.  This  is  obviously  very  important  because  edges  are  among  the  most 
important  features  we  use  to  interpret  images.  Blurring  is  also  the  most  limiting  property  of  many  noise 
reduction  techniques  especially  the  Fourier  based  methods.  Wavelet  based  noise  suppression  methods 
have  been  promising  largely  because  they  promise  to  leave  edges  unblurred.  The  spurious  peaks  at  the 
boundaries  between  monotonic  segments  are  the  only  deformation  of  the  signal  resulting  from  the  noise 
suppression.  It  is  caused  by  the  lack  of  averaging  at  the  ends  of  the  segments.  A  large  value  caused  by 
noise  in  the  center  of  the  segment  will  be  balanced  by  random  small  values  on  both  sides  of  it  while  a 
large  value  caused  by  noise  at  the  end  of  the  segment  satisfies  the  monotonic  criteria  and  is  not  changed. 
Another  way  to  see  the  effect  is  to  run  the  algorithm  over  a  series  of  random  values.  Some  part  of  the 
noise  energy  is  monotonic  and  it  is  primarily  at  the  ends  of  the  segment  as  shown  in  Figure  5.  That 
monotonic  part  of  the  random  noise  can  not  be  distinguished  from  real  monotonically  increasing  signal. 
Figure  2  shows  the  results  at  very  low  SNR’s.  The  noise  has  a  standard  deviation  of  one  and  the  signal 
peaks  have  a  height  of  one.  The  SNR  over  the  peaks  is  one  and  the  ratio  of  the  total  signal  energy  to  the 
total  noise  energy  is  0.6.  The  peaks  are  recovered  and  the  edges  are  sharp.  The  leakage  of  noise  at  the 
boundaries  between  the  segments  is  much  more  pronounced  than  at  higher  SNR,  as  might  be  expected. 
Figure  3  demonstrates  the  importance  of  extrema  selection.  All  three  peaks  were  recovered  when  the 
thresholds  were  properly  selected.  However,  if  the  threshold  for  extrema  selection  was  too  low,  too  many 
extrema  were  selected  and  extraneous  peaks  were  introduced  as  in  Figure  3B.  If  the  threshold  for  extrema 
selection  was  too  high,  too  few  extrema  were  selected  and  features  with  small  energy  were  lost  as  the  6 
pixel  wide  peak  was  lost  in  Figure  3D.  As  the  SNR  was  reduced  to  extremely  low  levels,  approximately 
0.1,  the  noise  leakage  at  the  boundaries  between  monotonic  segments  dominates  the  result  as  shown  in 
Figure  4.  Features  that  have  a  large  enough  area  were  still  found  and  the  sharp  edges  were  recovered 
effectively.  However,  the  noise  leakage  at  the  boundaries  between  monotonic  segments  overwhelms  the 
recovered  signal.  The  noise  leakage  at  the  boundaries  can  be  reduced  by  averaging  the  monotonic  fits 
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obtained  with  several  (5  to  11)  boundary  points  near  the  maximum.  If  the  average  is  weighted  by  the  error 
between  the  data  and  the  fit,  sharp  peaks  can  still  be  recovered.  However,  computation  times  increase 
enough  to  discourage  use  of  this  technique  on  images. 

III.  Methods 

The  first  method  of  processing  images  to  try  is  to  simply  process  the  rows  of  the  image,  then  rotate 
the  result  and  process  the  rows  again.  Then  iterate  by  rotating  the  result  and  processing  the  rows  until 
nothing  is  changing  and  quit.  This  simple  method  works  pretty  well  but  it  needs  to  be  modified  a  bit 
because  of  two  directional  effects. 

A.  Directional  Effects 

The  first  is  that  there  are  many  kinds  of  extrema  in  two  dimensions.  It  is  a  directional  property  in  two 
dimensions;  i.e.,  a  point  can  be  an  extrema  in  one  direction  but  not  in  another.  For  example,  a  saddle 
point  is  an  extrema  in  all  directions  but  two.  Therefore,  the  extrema  must  be  determined  independently 
for  each  direction  processed  in  an  image.  The  second  directional  effect  is  somewhat  more  subtle.  When  an 
image  is  processed  in  multiple  directions,  most  of  the  noise  is  suppressed  in  the  first  direction  processed. 
The  leakage  of  noise  at  the  boundary  between  monotonic  segments  is  greatest  perpendicular  to  the  first 
direction  processed  and  much  less  in  other  directions.  The  results  of  making  the  columns  monotonic  first 
and  then  the  rows  is  different  from  those  obtained  by  making  the  rows  monotonic  first  and  then  the 
columns.  If  this  effect  were  limited  to  the  edges  of  the  image  it  would  be  tolerable  but  it  occurs  at  each 
and  every  boundary  between  segments  throughout  the  image.  The  leakage  of  noise  at  boundaries  leads  to 
the  widening  of  a  peak  perpendicular  to  the  first  direction  made  monotonic.  The  effect  is  demonstrated 
in  Figure  6.  Just  noise  was  processed  to  show  how  it  is  blurred.  If  the  rows  are  made  monotonic  first,  the 
rows  stay  essentially  monotonic  after  the  columns  are  made  monotonic.  Therefore,  each  direction  needs 
to  be  made  monotonic  only  once. 

B.  Algorithm 

Our  method  of  processing  images  attempts  to  handle  these  directional  effects.  An  intermediate  image 
was  found  for  each  angle  used.  To  obtain  the  an  intermediate  image,  the  original  image  was  rotated 
to  the  appropriate  angle.  The  rotated  images  were  obtained  by  interpolating  to  a  rotated  grid  using 
cubic  B-splines  [9,  8].  The  rotated  version  of  the  image  was  blurred  using  a  simple  Gaussian  filter  in  the 
Fourier  domain.  The  extrema  along  the  rows  and  along  the  columns  of  the  blurred  image  were  identified. 
The  two  sets  of  extrema  were  examined  and  those  that  differed  from  adjacent  extrema  by  less  than  the 
threshold  were  eliminated.  Each  row  of  the  image  was  made  monotonic  between  extrema.  Then  the 
columns  of  the  row  processed  image  were  made  monotonic  between  extrema.  The  row-column  processed 
image  was  averaged  with  the  column-row  processed  image  to  obtain  the  intermediate  image.  Making 
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the  image  monotonic  in  two  directions  is  sufficient  to  suppress  most  of  the  noise  and  keep  most  of  the 
features.  Processing  the  image  at  more  angles  to  get  the  intermediate  image  takes  too  much  processing 
time.  Each  intermediate  image  had  the  noise  leakage  primarily  perpendicular  to  the  first  direction  made 
monotonic.  To  remove  the  noise  leakage,  intermediate  images  at  several  angles  were  averaged  together. 
Actually  the  intermediate  images  were  least  squares  fit  to  the  original  image  to  obtain  the  final  processed 
image.  A  simple  average  can  also  be  used  instead  of  a  least  squares  fit.  The  final  processed  image  is 
almost  monotonic  between  extrema  in  all  directions  and  contains  no  preferred  directional  effects. 

IV.  Results 

The  preliminary  results  are  very  promising.  Noise  can  be  suppressed  significantly.  Good  quality  images 
can  be  extracted  from  noise  that  would  otherwise  make  the  images  of  very  limited  usefulness.  The  first 
example  is  an  MRI  of  a  phantom  shown  in  Figure  7.  The  imaging  technique  selected  produced  very- 
noisy  images.  The  ‘true’  image  was  obtained  from  averaging  64  of  the  images  to  suppress  noise  through 
averaging.  Four  angles  were  used  in  the  noise  suppression.  The  edges  and  most  of  the  small  features 
are  recovered  well.  The  first  two  rows  of  small  pins  in  the  black  field  on  the  upper  left  section  of  the 
phantom  are  retained.  The  edges  remain  sharp.  However,  the  fan  shaped  resolution  pattern  is  less  well 
seen.  Smaller  parts  of  the  fan  pattern  could  probably  be  visualized  if  more  angles  were  used.  Most  of 
the  noise  remaining  in  the  noise  suppressed  image  is  leakage  at  the  boundaries  between  segments.  The 
second  example  is  of  a  functional  MRI  (fMRI)  image.  Functional  MRI  is  a  procedure  that  identifies  areas 
of  the  brain  that  are  used  to  perform  mental  tasks  such  as  the  motor  strip  or  the  areas  used  for  different 
kinds  of  memory.  Surgeons  need  such  information  to  avoid  critical  areas  of  the  brain  when  operating.  It 
is  also  important  to  understand  damage  and  reorganization  of  the  brain  following  injury  such  as  stroke 
or  trauma.  In  fMRI  many  sequential  images  of  the  brain  are  obtained  during  a  mental  task  and  during 
rest  from  that  task.  The  area  of  the  brain  that  is  active  during  the  task  has  very  small  signal  intensity 
changes  that  correlate  with  the  task  being  on  and  off.  These  signal  changes  are  caused  by  increased  blood 
flow  during  activation  and  are  only  around  two  percent  of  the  signal.  Noise  or  patient  movement  can 
render  the  study  inconclusive.  The  image. shown  in  Figure  8  is  a  single  image  from  an  fMRI  study.  The 
monotonic  noise  suppression  produces  excellent  results.  Eight  angles  were  used  in  the  noise  suppression. 
The  noise  suppressed  image  has  almost  all  the  features  present  in  the  original  but  the  noise  is  reduced 
significantly.  The  gray  and  white  matter  can  also  be  more  easily  segmented  in  the  noise  suppressed  image. 
The  regions  of  similar  signal  intensity  are  grouped  very  well.  The  third  example  is  an  T1  weighted  MR 
image  of  a  brain  study.  One  of  the  important  features  of  T1  weighted  images  is  to  differentiate  gray  and 
white  matter.  The  gray  matter  is  the  less  bright  shell  of  tissue  forming  the  cortical  surface.  The  brighter 
white  matter  is  inside  the  gray  and  the  dark  cerebrospinal  fluid  that  cushions  the  biain  is  between  the 
gray  matter  and  the  skull.  Eight  angles  were  used  in  the  noise  suppression.  Noise  was  added  to  the 
image  to  reduce  the  SNR  from  80  to  16  and  more  noise  was  added  to  decrease  the  SNR  farther  to  4.  The 
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monotonic  noise  suppression  increased  the  SNR  from  16  to  85  and  from  4  to  30.  The  gray-white  matter 
separation  is  recovered  well. 

A.  Applications 

There  are  many  potential  applications  of  monotonic  noise  suppression  but  two  applications  are  the  most 
promising.  Segmentation  is  the  most  natural.  Areas  that  have  similar  signal  intensities  are  grouped  into 
contiguous  regions  very  effectively.  Recovery  of  images  in  very  high  noise  is  the  other  natural  application 
because  this  technique  works  well  at  high  SNR’s  and  at  very  low  SNR’s.  Most  noise  reduction  methods  do 
not  perform  well  at  low  SNR’s.  For  example,  in  mammography  we  are  using  monotonic  noise  suppression 
as  a  preprocessing  step  for  a  watershed  algorithm  that  is  used  to  segment  the  digitized  mammograms. 
Monotonic  noise  suppression  can  be  adjusted  to  remove  the  small  features  that  tend  to  make  the  watershed 
algorithm  oversegment. 

jB.  Future  Improvements 

This  technique  is  still  being  actively  developed  and  there  are  many  improvements  to  be  made.  Many 
of  the  improvements  necessary  require  better  theoretical  understanding.  For  example,  how  many  angles 
need  to  be  used  to  suppress  noise  adequately?  The  answer  will  depend  on  the  area  of  the  regions  between 
extrema;  more  angles  will  be  needed  for  larger  areas.  How  fast  can  extrema  change  with  angle?  Most 
extrema  are  found  if  only  four  angles  are  used.  Can  the  leakage  at  the  boundaries  be  suppressed  without 
eliminating  smaller  features?  The  most  important  improvement  to  make  is  to  make  the  selection  of  extrema 
more  robust  and  automatic.  The  same  tricks  that  are  used  to  find  maximum  gradient  points  in  wavelet 
denoising  can  be  applied  here  as  well.  For  example,  tracking  extrema  across  multiple  scales  and  using 
continuity  would  probably  help.  Also  estimating  the  cutoffs  that  determine  the  extrema  automatically 
from  the  estimated  SNR  would  help.  Now  the  selected  extrema  are  shown  on  the  screen  to  select  the 
correct  cutoff  parameters. 

V.  Conclusions 

Forcing  monotonic  change  between  extrema  is  a  very  promising  noise  suppression  technique.  Edges  are 
not  blurred  and  it  works  well  over  a  wide  range  of  SNR’s.  As  the  SNR  drops,  features  that  have  less 
energy  than  the  noise  spikes  are  lost.  The  key  to  recovering  the  image  accurately  is  to  identify  the  correct 


extrema. 
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Fig.  1 .  An  example  of  suppressing  noise  on  a  line  segment  with  a  simulated  signal  and  normally  distributed 
random  noise.  The  signal  is  three  boxcars.  The  height  of  the  boxcars  is  1.  The  boxcars,  from  left  to 
right,  are  35,  45  and  6  pixels  wide.  The  signal  with  added  noise  is  shown  in  A.  The  standard  deviation 
of  the  normally  distributed  random  noise  is  0.1.  The  Gaussian  smoothed  signal  with  the  extrema 
selected  are  shown  in  B.  Standard  deviation  of  the  Gaussian  is  4  percent  of  the  total  bandwidth  and 
the  extrema  were  required  to  change  by  20  percent  of  the  maximum  signal.  The  result  of  forcing 
monotonic  change  between  the  extrema  is  shown  in  C.  The  original  signal  before  noise  was  added  is 
shown  in  D  for  comparison. 
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Fig.  3.  The  result  of  suboptimal  selection  of  extrema.  Figure  2  was  reconstructed  with  extrema  that 
were  larger  than  15  percent  of  the  maximum.  Figures  A  and  B  show  the  result  of  selecting  extrema 
larger  than  10  percent  of  the  maximum.  Figures  A  and  C  correspond  to  B  in  Figure  2  and  Figures 
B  and  D  correspond  to  C  in  Figure  2.  The  lower  threshold  result  in  two  new  extrema  seen  to  the  far 
left  of  the  segment.  That  extra  extrema  produces  a  rather  large  peak  in  the  final  result.  Figures  C 
and  D  show  the  result  of  selecting  the  extrema  larger  than  20  percent  of  the  maximum.  The  extrema 
over  the  6  pixel  wide  spike  in  the  signal  is  lost  and  does  not  appear  in  the  final  result. 


[ 


Fig.  4.  The  result  with  extremely  large  noise.  The  same  signal  as  in  Figure  1  was  corrupted  with  noise 
with  a  standard  deviation  of  10.  The  format  is  the  same  as  in  Figure  1.  Figures  A-D  are  scaled  the 


same. 


Fig.  5.  The  result  of  forcing  monotonic  change  on  normally  distributed  random  noise  shows  the  leakage 
of  noise  at  the  boundaries  between  monotonic  segments.  The  random  series  was  monotonically  fit 
to  two  segments;  the  first  increasing  and  the  second  decreasing.  There  is  always  a  false  peak  at  the 
boundaries  between  monotonic  segments  because  the  monotonic  fit  allows  negative  noise  to  pass  at 
the  beginning  of  the  segment  and  positive  noise  to  pass  at  the  end  of  an  increasing  segment.  In  the 
middle  of  the  segment,  values  to  either  side  are  enough  to  average  out  the  positive  and  negative  noise. 
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Fig.  6.  A  field  of  random  numbers  with  standard  deviation  1024  was  made  monotonic  along  the  columns 
first  and  then  along  the  rows.  The  center  row  were  the  maxima  when  the  columns  were  processed  and 
the  central  column  were  the  maxima  when  the  rows  were  processed.  The  result  shows  the  directional 
blurring  of  noise  perpendicular  to  the  direction  first  processed. 


Fig.  7.  Figure  A  shows  an  MRI  image  of  a  phantom.  The  imaging  technique  used  produced  a  very  poor 
SNR  image,  shown  in  Figure  A.  Figure  B  shows  the  results  of  noise  reduction  using  monotonic  fits 
between  extrema.  Figure  C  shows  a  Fourier  filtered  image  with  the  same  SNR  as  in  Figure  B;  the 
blurring  limits  the  usefulness  of  the  image.  Figure  D  shows  the  average  of  64  acquisitions  averaged 
together  to  reduce  noise.  The  SNR  in  Figure  D  is  8  times  that  in  Figure  A. 


Fig.  8.  Figure  A  shows  the  original  axial  MR  image  and  Figure  B  shows  the  results  of  forced  mono¬ 
tonic  noise  suppression.  The  parameters  were  set  to  produce  uniform  areas  that  could  be  useful  in 
segmentation.  The  SNR  is  improved  significantly  by  forced  monotonic  noise  suppression. 


Fig.  9.  Figure  A  shows  the  original  T1  weighted  sagital  MR  image  that  has  a  relatively  high  SNR, 
SNR=80  within  the  brain.  Figures  B  and  C  show  the  original  image  with  added  normally  distributed 
random  noise,  SNR=:16  and  SNR=4  within  the  brain  for  Figures  B  and  C  respectively.  Figure  D 
shows  the  Gaussian  filtered  image  with  the  selected  extrema  superimposed.  Figure  E  shows  the  result 
of  forced  monotonic  noise  suppression  on  Figure  B.  The  SNR  increased  from  16  to  85.  Figure  F  shows 
the  result  of  monotonic  noise  suppression  on  Figure  C.  The  SNR  is  increased  to  from  4  to  30. 


Monotonic  Noise  Suppression  Used  to  Improve 
the  Sensitivity  of  fMRI  Activation  Maps 

John  B.  Weaver 


We  have  introduced  a  new  method  of  removing  noise 
from  images  that  identifies  significant  extrema  and 
forces  the  pixel  intensities  between  any  two  extrema 
to  change  monotonicly.  The  method  has  some  similari¬ 
ties  to  wavelet  denoising  methods  we  worked  with 
several  years  ago  but  is  generally  more  stable  and  is 
effective  on  images  with  lower  SIMR's.  In  this  paper  the 
method  of  monotonic  filtering  is  used  to  increase  the 
sensitivity  in  functional  magnetic  resonance  imaging 
(fMRI)  studies.  We  have  used  the  increased  sensitivity 
to  improve  the  temporal  resolution  in  fMRI  studies  by 
roughly  a  factor  of  six.  A  motor  activation  study  was 
acquired  with  single  slice  256  x  256  pixel  T2*  weighted 
images;  six  cycles  of  finger  tapping  were  acquired. 
Each  cycle  consisted  of  five  images  of  rest  followed  by 
five  images  of  right  hand  finger  tapping  followed  by 
five  images  of  left  hand  finger  tapping.  The  z-scores 
were  calculated  and  used  as  the  activation  map.  The 
left  and  right  activations  were  both  clearly  visible 
when  all  six  cycles  were  used  in  the  analysis.  However, 
no  definitive  activation  was  seen  for  any  one  cycle. 
When  the  original  256  x  256  images  were  averaged 
down  to  64  x  64  pixel  images  before  calculation  of  the 
z-scores,  the  activations  were  partially  identified.  When 
the  original  images  were  filtered  using  the  monotonic 
noise  reduction  algorithm,  the  left  activation  was 
clearly  visible  in  three  of  the  six  cycles  and  partially 
visible  in  two  others.  The  right  activation  was  partially 
visible  in  4  out  of  6  cycles.  Optimized  noise  reduction 
should  improve  the  results  significantly.  The  ability  to 
use  a  single  cycles  is  very  important  in  fMRI  studies 
because  many  stimuli  are  more  difficult  to  maintain 
over  many  cycles  and  because  complex  processes 
such  as  in  cognitive  or  memory  activity  do  not  have 
simple  responses. 

Copyright©  1998  by  W.B.  Saunders  Company 

Noise  reduction  has  many  applications: 

eg,  as  a  preprocessing  step  for  many  image 
interpretation  algorithms,  in  segmentation,  in  con¬ 
trast  enhancement  among  many  others.  Noise  reduc¬ 
tion  can  be  seen  as  finding  the  features  in  the  image 
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and  preferentially  separating  the  features  from  the 
noise  in  some  way.  Fourier  filters  reduce  the  high 
frequencies  where  noise  is  assumed  to  be  dominant 
and  leave  the  low  frequencies  where  signal  is 
assumed  to  be  dominant.  The  undesired  result  is 
bandlimited  blurring  of  real  features.  Wavelet  de- 
noising‘‘^  improves  the  situation  by  allowing  the 
bandwidth  passed  to  be  increased  around  edges  in 
the  image  and  decreased  in  flat  regions.  This 
mitigates  the  bandlimited  blurring  if  the  edges  can 
be  found  in  the  image.  Indeed  the  key  to  wavelet 
denoising  is  correctly  identifying  the  edges.  Shrink¬ 
age  algorithms  assume  that  large  wavelet  coeffi¬ 
cients  are  likely  to  be  a  result  of  an  edge  and 
smaller  wavelet  coefficients  are  likely  to  be  noise. 
Therefore,  the  larger  coefficients  are  kept  and  the 
smaller  ones  are  suppressed. 

Functional  magnetic  resonance  imaging  (fMRI) 
looks  for  small  changes  in  blood  flow  that  are 
produced  by  mental  activity  associated  with  mental 
tasks.  The  signal  changes  are  on  the  order  of  the 
noise  in  MR  images.  They  are  generally  around 
3%P  Identifying  these  changes  in  MR  images  with 
comparable  noise  is  a  difficult  task.  The  signal  can 
only  be  found  by  averaging  over  many  cycles  of  the 
task.  However,  many  tasks  are  difficult  to  consis¬ 
tently  sustain  over  multiple  cycles;  e.g.,  cognitive 
tasks  or  other  tasks  with  multiple  parts  and  com¬ 
plex  activations.  We  are  using  monotonic  noise 
reduction  to  reduce  the  number  of  activation  cycles 
required  to  identify  the  region  activated. 

METHODS 

Monotonic  Noise  Reduction 

Monotonic  noise  rediiction^-^  finds  the  significant  extrema, 
both  the  maxima  and  the  minima,  and  forces  monotonic  change 
between  the  extrema.  The  procedure  is  best  described  in  one 
dimension  first  and  then  extended  to  two  dimensions.  The 
significant  features  are  the  extrema  in  a  slightly  blurred  image 
that  are  different  from  adjacent  extrema  by  a  user  defined 
minimum  value.  The  .selection  of  the  correct  extrema  is  of 
critical  importance.  W^hen  the  significant  extrema  have  been 
selected,  the  signal  must  by  definition  change  monotonically 
between  any  pair  of  adjacent  extrema.  Monotonic  change  is 
forced  by  using  an  elegant,  simple  algorithm  for  finding  the 
monotonic  series  of  numbers  that  fits  the  original  series  best  in 
the  least  squares  sense. The  monotonic  fit  averages  large 
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Figl.  A  flow  chart  of  the  steps  taken 
to  obtain  the  extrema  in  ail  four  direc¬ 
tions  that  are  used  in  the  monotonic 
noise  reduction  algorithm.  The  image  is 
blurred  using  a  simple  Gaussian  filter  In 
the  Fourier  domain  to  average  out  small 
extrema.  Then  the  pixels  are  reordered 
into  strings  that  follow  the  four  primary 
directions  and  the  extrema  are  selected 
In  each  reordered  string  of  numbers. 
Only  extrema  that  are  different  from 
adjacent  extrema  by  more  than  a  thresh¬ 
old  are  included  in  the  four  final  sets  of 
extrema.  The  user  selects  the  width  of 
the  Fourier  filter  and  minimum  differ¬ 
ence  between  extrema.  These  two  pa¬ 
rameters  determine  how  small  a  feature 
Is  rejected  as  noise  or  is  Included  in  the 
final  image. 


Fig  2.  A  flow  chart  of  the  steps  taken  in  monotonic  noise  reduction.  The  pixels  in  the  original  image  are  reordered  along  the  four 
primary  directions  and  the  pixel  values  between  extrema  are  fit  to  monotonic  functions.  The  resulting  series  of  pixel  values  are 
reordered  in  the  direction  orthogonal  to  the  first  direction  and  the  pixel  values  between  that  set  of  extrema  are  fit  to  monotonic 
functions.  The  resulting  four  images  are  averaged  to  form  the  final  filtered  image. 
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Fig  3.  The  fifth  resting  image  in  the  raw  data  from  a  motor 
activation  study. 


Fig  5.  The  fifth  image  in  the  motor  study  after  It  has  been 
averaged  down  to  a  64  x  64  image.  Each  pixel  is  the  average 
of  sixteen  pixels  In  the  full  resolution  image.  The  noise  is 
reduced  but  the  resolution  suffers. 


enough  group.s  of  adjacent  points  in  the  series  to  achieve  a 
monotonic  progre.ssion.  It  is  robust  and  relatively  fast.  The 
resulting  monotonic  .series  is  not  ju.st  a  fit  to  a  limited  .set  of  basis 
functions.  The  monotonic  series  is  not  bandlimited  .so  edges  are 
not  blurred. 

The  extension  to  two  dimensions  is  complicated  because  an 
extrema  in  one  direction  is  not  necessarily  an  extrema  in  another 
direction.  For  example,  pixels  along  a  ridge  will  be  extrema 
perpendicular  to  the  ridge  but  not  along  the  ridge.  The  method 
we  use  is  to  find  the  extrema  in  four  directions  as  in  Fig  I  and 
filter  the  image  separately  in  four  directions:  horizontal,  vertical, 
diagonal  upper-left  to  lower-right,  and  diagonal  lower-left  to 
upper-right.  There  are  four  intermediate  images  formed  by  first 
filtering  the  original  image  in  each  of  the  four  primary  direc¬ 
tions.  The  intermediate  images  are  then  filtered  in  the  direction 
perpendicular  to  the  first  filter  direction.  Filtering  in  other 
directions  can  be  done  but  produces  little  change.  Once  an 
image  has  been  filtered  in  two  orthogonal  directions,  it  is 
essentially  monotonic  in  the  other  directions,  too.  However,  the 
four  intermediate  images  are  all  different  valid  solutions  that  are 
essentially  monoionic  in  all  four  directions.  We  average  the  four 


Fig  4.  The  results  of  monotonic  noise  reduction  on  the  fifth 
image  in  the  motor  activation  study.  The  noise  Is  significantly 
reduced  from  that  in  Fig  3. 


to  obtain  the  final  filtered  image.  The  procedure  is  outlined  in 
Fig  2. 

We  have  tried  a  weighted  least  squares  fit  of  tlie  four  to  the 
original  image  but  it  is  a  little  more  unstable  and  does  not 
produce  a  noticeably  belter  image.  We  have  also  filtered  the 
image  in  more  than  four  directions  but  the  most  improvement  in 
image  quality  is  obtained  in  the  first  four  directions  so  we 
usually  limit  ourselves  to  four  directions. 

fMRl  Protocol 

An  axial  plane  four  cm  below  the  apex  of  the  cranium  was 
imaged.  Six  cycles  of  rest  and  motor  activation  were  obtained. 
Each  cycle  of  images  consisted  of  five  resting  images  followed 
by  five  images  taken  during  tapping  of  the  right  hand  followed 
by  five  images  taken  during  tapping  of  the  left  hand.  Each 
gradient  echo  image  was  a  256  X  256,  24  cm  field  of  view, 
image  with  TR  of  70  ms  and  TE  of  40  ms.  The  images  were 
strongly  T2=‘'  weighted  to  identify  BOLD  contrast  changes. 


Fig  6.  The  average  of  all  90  images  from  the  study  provides 
a  relatively  noise  free  comparison.  The  noise  and  resolution 
resemble  that  in  the  monotonically  filtered  image  in  Fig.  4 
more  closely  than  that  in  the  original  Image  in  Fig  3  or  that  in 
the  reduced  resolution  image  In  Fig  5. 
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Fig  7.  The  full  resolution  256  x  256  pixel  reference  activation  maps  for  the  right  hand  motor  task  (on  the  left)  and  the  left  hand 
motor  activation  task  (on  the  right).  The  map  is  the  correlation  with  the  right  hand  motion  over  the  standard  deviation  of  the 
intensities  at  that  pixel  position.  All  six  cycles  were  used  to  obtain  the  activation  maps.  The  activations  are  clearly  present  in  the 
expected  locations. 


Activation  Maps 

Simple  z-.scores  for  each  pixel  were  used  to  identified 
activated  regions.  The  standard  deviation  of  the  ninety  images 
was  used  as  a  measure  of  the  total  variation  in  signal.  To  obtain 
the  reference  activation  maps  all  six  cycles  were  used  to  obtain 
the  z-scores.  To  obtain  the  signal  coirelated  with  right  hand 
finger  tapping,  the  average  of  the  rest  and  left  hand  tapping  was 
subtracted  from  the  average  of  the  right  hand  finger  tapping.  The 
activation  map  for  the  right  hand  motor  task  was  the  ratio  of  the 
signal  correlated  with  right  hand  finger  lapping  over  the  total 
variation  in  signal.  Similarly,  the  activation  map  for  the  left  hand 
motor  task  was  the  ratio  of  the  signal  correlated  with  left  hand 
finger  tapping  over  the  total  variation  in  signal.  No  realignment 
of  the  images  was  done. 


Single  Cycle  Tests 

The  reference  activation  maps  were  obtained  with  all  six  task 
cycles.  The  reference  activation  maps  were  used  as  the  gold 
standard  showing  the  true  regions  activated  which  were  taken  to 
be  the  motor  centers.  During  any  single  task  cycle  there  could  be 
transient  activation  in  other  areas  or  lack  of  activation  in  the 
regions  identified  in  the  reference  activation  maps.  However,  the 
regions  shown  in  the  reference  activation  maps  do  show  the 
mo.st  likely  areas  activated. 

The  activation  maps  for  each  cycle  were  calculated  and 
compared  to  the  reference  activation  maps.  The  images  were 
processed  in  two  ways  to  average  out  noise  in  the  images  and 
help  make  single  cycle  activations  definitive.  First,  the  images 
were  averaged  dou  n  from  256  by  256  pixels  to  64  by  64  pixels. 


Fig  8  The  full  resolution  256  x  256  pixel  activation  map  obtained  using  the  fourth  task  cycle.  The  right  hand  motor  task  is  on  the 
leh  and  the  left  hand  motor  activation  task  is  on  the  right.  No  activation  can  be  identified  without  reference  to  the  maps  in  Fig  7, 
Lere  are  a  few  pixels  in  the  motor  areas  but  not  enough  to  differentiate  the  motor  centers  from  other  areas  in  the  maps. 
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Fig  9.  The  ectivation  maps  obtained  from  the  full  resolution,  monoton.cally  f.ltered  .mages 
used  The  right  hand  motor  task  is  on  the  left  and  the  left  hand  motor  activation  task  is  on  the  right.  The  left  hand  '*  *= 

Vis  ble  but  more  dispersed  than  in  the  reference  activation  maps.  This  could  be  the  real  activation  »'’^™;"XaVb"t  h  Is  ^ot 
composite  activation  maps  or  it  could  be  the  result  of  noise.  The  right  hand  activation  is  visible  when  compared  to  Fig  7  but  ,t  ,s  not 

as  clear  as  the  left  hand  activation. 


The  averusing  reduced  noise  but  also  reduced  the  resolution. 
Second,  the  images  were  processed  with  monotonic  noise 
reduction.  The  resolution  remained  the  same  as  in  the  original 
images. 

RESULTS  AND  DISCUSSION 

The  fifth  image  of  the  study  is  shown  in  Fig  3. 
The  same  image  after  monotonic  noise  reduction 
and  after  averaging  are  shown  in  Figs  4-5.  The 
average  image  of  the  90  image  study  is  shown  in 


Fig  6  for  comparison.  The  noise  that  is  apparent  in 
Fig  3  is  much  reduced  in  Fig  4  without  loss  of 
significant  features.  The  resolution  in  Fig  5  suffers 
significantly. 

Figure  7  shows  the  reference  activation  maps 
obtained  with  all  six  task  cycles.  The  right  and  left 
motor  centers  are  clearly  visible.  The  locations  of 
both  motor  centers  are  also  in  roughly  the  same 
locations  as  in  most  subjects.  The  signal  that 
correlates  with  the  activation  is  around  seven  times 


visible.  The  right  hand  activation  is  not  visibly  different  from  other  areas  in  the  map. 
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Fiq  1 1 .  The  full  resolution  256  x  256  pixel  activation  map  obtained  using  the  fifth  task  cycle.  The  right  hand  motor  task  is  on  the 
left  and  the  left  hand  motor  activation  task  is  on  the  right.  No  activation  can  be  identified.  Regions  of  activation  can  only  be  seen  if 
compared  to  the  known  activation  maps  in  Fig  7. 


the  standard  deviation  in  the  activated  pixels  so  the 
confidence  in  the  activations  is  high. 

The  fourth  and  fifth  single  cycle  activation  maps 
are  shown  in  Figs  8-13.  The  other  single  cycle 
maps  are  similar.  No  definitive  activation  was  seen 
for  any  one  cycle  without  filtering  or  averaging  the 
images  prior  to  calculation  of  the  activation  map. 
The  activations  in  Figs  8  and  11  can  be  seen  in 
retrospect  if  compared  to  the  six  cycle  activations 
in  Fig  7.  The  activation  maps  obtained  from  the 
monotonically  filtered  images  in  Figs  9  and  12 
show  the  left  hand  activation  very  clearly.  The  right 


hand  activation  can  be  seen  retrospectively  if 
compared  to  the  six  cycle  activation  maps  in  Fig  7. 
The  reduced  resolution  activation  maps  in  Figs  10 
and  13  are  similar  to  the  monotonically  filtered 
activation  maps.  The  left  hand  activations  are 
visible  but  the  activation  for  the  right  hand  task  is 
even  less  visible  than  in  the  maps  obtained  from  the 
monotonically  filtered  images.  The  reduced  right 
hand  visibility  might  be  a  partial  volume  effect. 

A  method  to  optimize  the  monotonic  filter  set- 
tinsjs  is  needed  to  improve  the  results.  The  filter 
settings  were  set  arbitrarily  to  what  produced  a 


Fia  12  The  activation  maps  obtained  from  the  full  resolution,  monotonically  filtered  images  when  only  the  fifth  task  cycle  was 
used  The  right  hand  motor  task  is  on  the  left  and  the  left  hand  motor  activation  task  is  on  the  right.  The  left  hand 
visible  but  again  it  is  more  dispersed  than  in  the  reference  activation  maps.  The  right  hand  activation  is  also  only  visible  when 

compared  to  Fig  7. 
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Fig  13,  The  64  by  64  pixel,  activation  maps  obtained  from  the  reduced  resolution  images  when  only  the  fifth  task  cycle  was  used. 
The  right  hand  motor  task  is  on  the  left  and  the  left  hand  motor  activation  task  is  on  the  right.  The  left  hand  activation  is  clearly 
visible.  The  right  hand  activation  is  not  visibly  different  from  other  areas  in  the  map. 


visually  good  image.  Better  results  can  probably  be 
obtained  with  optimal  filter  settings. 

CONCLUSIONS 

Forcing  monotonic  change  between  extrema  is  a 
promising  noise  suppression  technique.  Edges  are 
not  blurred  and  it  works  well  over  a  wide  range  of 
SNR's.  As  the  SNR  drops,  features  that  have  less 
energy  than  the  noise  spikes  are  lost  but  edges  are 
not  bluiTed.  The  key  to  recovering  the  image 
accurately  is  to  identify  the  correct  extrema. 

We  used  the  monotonic  filtering  method  to 
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reduce  the  noise  in  functional  magnetic  resonance 
(fMRI)  images.  It  improves  the  SNR  enough  to 
identify  motor  activations  from  a  single  task  cycle 
rather  than  multiple  cycles  in  some  cases.  The 
increased  visibility  of  activated  regions  is  compa¬ 
rable  to  that  obtained  by  averaging  the  pixels  in  the 
images  to  reduce  the  noise.  However,  the  resolution 
is  not  degraded  when  monotonic  noise  reduction  is 
used  as  it  is  when  the  pixels  are  averaged.  Because 
monotonic  filtering  is  effective  at  much  lower 
SNR's  than  the  common  wavelet  filtering  methods, 
it  is  particularly  useful  in  fMRI  studies  where  the 
noise  is  relatively  high. 
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Tissue  elasticity  is  thought  to  hold  great  promise  for  ‘ 
the  diagnosis  of  cancer  because  It  is  often  harder  than 
normal  tissue.  MR  based  methods  encode  the  static 
displacements  [1]  or  dynamic  displacements  [2,3]  in  the 
phase  and  reconstruct  the  elasticity  from  the  motion. 

We  have  used  dynamic  displacements  because 
density  and  frequency  effects  can  be  evaluated  as  well 
as  the  elasticity.  The  gradient  echo  sequence  below  was 
implemented  on  a  GE  MRI.  The  KF  amplifier  was 
blanked  during  the  lOOHz  sinusoid  on  the  RF  channel. 
It  was  fed  into  the  power  amplifier  for  the  piezoelectric 
actuator  stacks.  The  lOOHz  sinusoidal  gradients 
encoded  the  harmonic  displacement.  The  3  gradients ' 
allowed  displacement  in  all  3  directions  to  be  measured. 
The  relative  phase  between  the  gradients  and  the 
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^^otion  Encoding  (ir^ients 
motion  was  controlled  very  precisely  by  changing  the'; 
starting  time  of  the  gradient  pulse.  Four  starting  times ' 
were  used  and  phase  at  each  position  was  fit  to: 


0  =  C  +  A  cos(ct)t  +  P) .  The  first  term  is  an  irrelevant 
constant.  Tlie  amplitude  and  phase  of  the  cos  term 
describe  the  harmonic  motion  completely. 

The  motion  was  produced  by  three  stacks  of 
piezoelectric  actuators..  Three  actuators  were  used  to 
get  sufficient  displacement.  Three  stacks  gave 
mechanical  stability.  The  system  was  current  limited 
with  one  audio  amplifier;  in  Independent  measurements 
the  device  produced  24  micron  displacements.  The 
specifications  state  that  the  stacks  can  provide  42 
micron  displacements  with  sufficient  current.  The 
■device  generates . large  forces  which  is  an  important 
■  advantage. 


Discussion:  Because  the  measured  displacement 
patterns  match  the  calculated  patterns  quite  well,  we 
are  confident  that  we  have  made  accurate 
measurements  of  the  motion  and  accurate  calculations 
of  the  motion.  The  magnitude  of  the  simulated 
displacements  are  larger  than  the  measured  values. 
Damping,  transient  effects,  and  changes  in  Poisson’s . 
ratio  could  all  account  for  the  variation. 
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and  Y  directions  in  a  homogeneous  agar  gel  are  shown 
in  Figs,  a  and  b.  The  TR=300ms.  TE=58ms.  NEX=2, 
■  SL=lcm,  FOV=8cm,  resolution  is  128  x  256.  Figs,  c 
and  d  show  the  forward  simulation  for  that  geometry 
.  and  material. 
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Abstract 

Recently,  imaging  modalites  such  as  ultrasound  and  magentic 
resonance  (MR)  have  been  used  to  measure  subsurface  dis¬ 
placements  in  tissue  and  several  inversion  schemes  have  been 
proposed  to  solve  for  stiffness  properties.  We  have  developed 
a  finite  element  based  inversion  scheme  which  operates  in  a 
sweeping  fashion  on  small  overlapping  subzones  of  the  tissue 
space.  The  zone  approach  allows  for  a  high  degree  of  spatial 
discretization  while  maintaining  algorithm  convergence.  Addi¬ 
tionally,  we  are  using  a  harmonic  elastodynamic  tissue  model 
as  the  basis  of  our  inversion  and  have  shown  accurate  recon¬ 
struction  simulations  with  up  to  15%  added  noise. 

Introduction 

Palpation,  although  effective  at  diagnosing  large  near-surface 
cancerous  tissue,  is  not  an  adequate  technique  for  detecting 
small  deep  tumors.  However,  based  on  the  success  of  palpa¬ 
tion,  the  high  contrast  in  stiffness  between  healthy  and  can¬ 
cerous  tissues  remains  an  impetus  for  developing  an  elasto- 
graphic  imaging  modality.  MR  and  ultrasound  elastography 
are  the  first  imaging  modalities  to  provide  subsurface  displace-, 
ment  data,  which  consequently  provides  tissue  strain  informa¬ 
tion,  that  can  be  used  in  an  inverse  method  to  recover  stiffness 
properties. 

Methods 

Previous  work  has  focused  on  estimating  stiffness  properties  by 
calculating  local  wavelengths  resulting  from  shear  excitation  of 
the  tissue  [1].  More  recently  model  based  reconstruction  has 
been  used  in  this  same  context  [2].  It  has  been  our  experience 
that  shear  wave  excitation  of  tissue  is  limited  due  to  attenu¬ 
ation.  In  our  approach,  we  use  longitudinal  harmonic  waves 
which  ultimately  generate  standing  shear  waves  deep  within 
the  tissue. 

The  equations  describing  the  elastodynamic  response  of  soft 
tissue  under  an  applied  harmonic  deformation  are, 


where  r  is  the  stress  tensor,  u  is  the  displacement  vector  and 
p  is  the  tissue  density.  Assuming  that  the  material  is  excited 
harmonically  at  frequency  cj,  these  equations  can  be  solved  in 
the  frequency  domain, 

=  V  •  T.  (2) 

where  u  =  Ue’"*.  Depending  on  th^  constitutive  relationships 
assumed  between  stress  and  strain,  the  presence  of  damping 
can  be  incorporated;  however,  for  this  discussion  we  have  as¬ 
sumed  Hookean  dependence  with  a  constant  Poisson’s  ratio 
of  1/  =  0.49;  thus  leaving  Young’s  modulus  the  only  unknown 
in  the  domain  (recall  that  displacement  data,  u,  is  measured 
from  the  imaging  technique). 

The  inversion  problem  is  a  nonlinear  Newton-based  itera¬ 
tive  scheme  which  minimizes  the  square  of  the  error  between 
measured  and  model-predicted  values  for  each  zone  and  solves 
for  the  distribution  of  Young’s  modulus.  The  zone  domain  is 
radially  shaped  and  determined  by  a  hierarchial  ordering  of 
local  residual  errors  which  cover  the  entire  mesh  (zone  bound¬ 
ary  conditions  are  determined*  from  the  MR  dataset).  After 


all  areas  of  the  mesh  have  been  iterated  on  a  specified  number 
of  times  (zone  iterations  vary  due  to  overlapping  capability),  a 
global  forward  problem  is  executed  using  the  updated  modulus 
distribution  and  the  zone  process  begins  again.  The  main  ad¬ 
vantage  of  this  technique  is  that  it  allows  sufficient  discretiza¬ 
tion  to  resolve  the  wavelengths  found  in  soft  tissue  harmonic 
motion. 

Results  &:  Discussion 

Using  this  inversion  scheme,  a  simulation  was  performed  on 
a  breast  cross-section  with  complex  phantom  shapes  of  vary¬ 
ing  contrast. '  Random  noise  (up  to  15  %  of  original  displace¬ 
ment)  was  added  to  the  forward  solution  data  to  simulate  sig¬ 
nal  degradation  in  the  MR  measurements.  This  noisy  solution 
is  then  operated  on  by  the  inversion  algorithm  described  pre¬ 
viously.  Figure  1  shows  the  results  with  an  initial  guess  of  a 
uniform  Young’s  modulus  of  7000  Pa.  The  inversion  process 


b 

Figure  1:  Breast  computational  phantom/reconstruction  with 
regions  of  varying  contrast  (2x,  5x,  lOx):  (a)  phantom  Young’s 
modulus  distribution,  (b)  inverse  solution  with  15%  random 
noise  added  to  data. 

consisted  of  18  sweeps  over  the  entire  space,  each  sweep  using 
roughly  1000  zones  of  about  150  elements  and  100  nodes  to 
insure  that  every  node  within  the  discretization  was  operated 
on  at  least  once.  Overall  the  results  shown  in  our  simulations 
are  extremely  encouraging.  This  work  was  supported  by  NIH 
grant  R01-NS33900  awarded  by  the  NINDS. 

References 

1.  R.  Mathupillai,  et  al.  Magnetic  Resonance  in  Medicine,  36, 
266-274,  1996. 

2.  A.  Manduca,  et  al.  Lecture  Notes  in  Computer  Science, 
1496,  606-613,  1998. 


AAPM  1998 


Medical  Physics  25(7)  Part  1,  P.  A 122-3 
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We  introduced  monotonic  noise  suppression  and  applied  it  to  improving  fMRI  activation 
maps.  Now  we  are  using  monotonic  noise  suppression  to  increasing  contrast  in  images.  Almost  any 
noise  reduction  algorithm  can  be  used  to  enhance  contrast.  Monotonic  noise  suppression  has  some 
advantages  over  the  wavelet  denoising  based  contrast  enhancement  method  we  introduced 
previously. 

The  basic  idea  is  to  subtract  two  versions  of  the  image  obtained  by  using  two  different  noise 
thresholds.  The  first  image  is  filtered  to  remove  just  the  noise.  The  second  image  is  filtered  to 
remove  noise  and  small  features  leaving  only  the  large  features.  When  the  second  image  is 
subtracted  from  the  first,  the  large  features  are  subtracted  away  leaving  the  small  features.  Actually 
the  second  image  times  some  factor  less  than  one  is  subtracted;  so  the  large  features  are  reduced 
rather  than  eliminated.  The  amount  of  contrast  enhancement  can  be  changed  by  changing  the  linear 
combination  of  filtered  images. 

Two  features  of  the  monotonic  noise  reduction  method  make  it  attractive  for  this  application. 
First,  it  is  not  band-limited.  Edges  are  not  blurred  at  all.  There  is  no  reduction  of  sharpness  in  the 
contrast  enhanced  image.  Second,  there  is  no  ringing  caused  by  undersampling.  When  a  noise 
reduction  method  that  introduces  ringing  into  the  image,  the  ringing  can  be  amplified  and  distorted 
to  look  like  small  features  in  the  contrast  enhanced  image.  Fourier  noise  reduction  methods  are 
infamous  for  introducing  ringing  and  wavelet  denoising  also  introduces  ringing  although  to  a  lesser 
extent. 
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Contrast  enhancement  requires  the  balance  of  two  competing  features:  making  small  changes  in 
the  image  bigger  and  controlling  the  growth  of  noise.  There  are  two  general  types  of  techniques  that 
have  been  used  in  the  past  to  enhance  contrast.  The  first  is  contrast  equalization.  However,  it  does 
not  do  well  when  there  are  big  bright  and  big  dark  areas  with  low  contrast  objects  in  each  because 
the  large  bright  and  dark  areas  dominate  the  histogram  and  do  not  help  the  low  contrast  objects  in 
between.  (See  Gonzalez  and  Wintz,  Digital  Image  Processing,  Section  4.2.3) 

The  second  class  of  contrast  enhancement  techniques  that  have  been  used  in  medical  imaging  are 
variants  of  transform  based  noise  reduction  methods.  We  developed  a  wavelet  based  algorithm  that 
is  used  in  mammography  by  several  groups  [Lu  1994].  One  way  of  seeing  the  process  is  that  the 
coarse  shape  of  the  object  (formed  by  the  coarse  scale  coefficients)  is  extracted  and  allowed  less  of 
the  dynamic  range  to  make  room  for  the  midscale  coefficients.  The  idea  is  the  same  as  the  old 
blurred  mask  or  unsharp  masking  methods. 

The  contrast  enhancement  method  we  are  presenting  here  follows  the  same  form  as  the  wavelet 
based  one  described  above.  But,  we  are  using  the  monotonic  noise  reduction  algorithm  [Weaver 
1997  a&b]  that  we  have  used  for  noise  reduction  in  fMRI  studies  [Weaver  1998]. 

We  describe  monotonic  noise  reduction  and  then  shown  how  it  is  used  to  enhance  contrast. 

Monotonic  Noise  Reduction: 


n 


Figure  :  Monotonic  noise  reduction  works  in  ID.  The  layout  is  the  same  for  both  figures.  The  top  is  the  original  signal. 
The  next  line  is  the  signal  with  added  noise;  SNR=10  on  the  left  and  SNR=1  on  the  right.  The  next  line  is  the  blurred 
signal  with  the  selected  extrema.  The  filtered  signal  is  shown  at  the  bottom.  The  filtered  signal  is  the  best  monotonic  fit 
between  selected  extrema  The  results  are  quite  good  even  at  very  low  SNR.  There  is  no  blurring  of  the  edges  and  no 
rineine.  The  only  distortion  is  the  noise  leakage  at  the  extrema. 


Figure  :  A  sagital  MRI  image  is  shown  at  the  left;  SNR=80.  The  middle  is  the  original  image  with  added  noise: 
SNR=16.  The  right  is  the  result  of  monotonic  noise  reduction;  SNR=85.  The  images  are  of  good  quality  with  no 
blurring  or  distortion.  Monotonic  noise  reduction  works  with  much  lower  SNR  that  common  wavelet  denoising 
algorithms. 


AAPM  1998 


Medical  Physics  25(7)  Part  1,  P.  A 122-3 


Contrast  Enhancement:  Two  features  of  the  monotonic  noise  reduction  method  make  it 

attractive  for  contrast  enhancement.  1)  It  is  not  band-limited;  edges  are  not  blurred  at  all.  So  there  is 
no  reduction  of  sharpness  in  the  contrast  enhanced  image.  2)  There  is  no  ringing  caused  by 
undersampling  in  monotonic  noise  reduction.  Ringing  must  be  controlled  just  as  noise  must  in  the 
contrast  enhanced  images  or  they  will  dominate  the  image. 


Figure  1 :  An  image  from  an  fMRI  study.  The  contrast  is  low  as  it  generally  is  in  fMRI  studies.  The  image  on  the  right 
is  the  contrast  enhanced  image.  The  darker  areas  are  darker  and  the  brighter  areas  are  brighter  than  in  the  original  image. 
No  visible  artifacts  were  introduced  in  the  contrast  enhanced  image.  However,  the  artifacts  that  were  present  in  the 
original  image  are  enhanced  just  as  the  small  features  are.  The  general  intensity  distribution  remains  the  same  so  the 
"look  and  feel"  of  the  contrast  enhanced  image  is  the  same  as  the  original  image. 


Figure  2:  The  top  curve  is  the  180th  column  of  the  original  image.  The  bottom  curve  is  the  180th  column  of  the  contrast 
enhanced  image.  The  same  peaks  are  present  but  the  size  of  the  peaks  are  equalized:  the  larger  discontinuities  are 
relatively  smaller  and  the  smaller  peaks  are  relatively  larger.  The  relative  sizes  of  the  peaks  are  generally  the  same  so  the 
"look  and  feel”  of  the  image  is  the  same  unlike  histogram  equalized  images. 

Figs.  1  &  2  show  that  the  contrast  has  been  enhanced  without  allowing  noise  to  dominate.  Each 
feature  in  the  contrast  enhanced  image  can  be  seen  in  the  original  image  but  less  clearly. 
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Elasticity  estimates  using  phase  contrast  MRI  measurements  of 
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by 

John  B.  Weaver  f,  Elijah  van  Houten  §,  Michael  I.  Miga  §,  Francis  E.  Kennedy  §,  Steven 
P.  Poplack  t,  Helene  M.  Nagy  f,  Keith  D.  Paulsen  § 

t  Department  of  Radiology,  Dartmouth-Hitchcock  Medical  Center 
§  Thayer  School  of  Engineering,  Dartmouth  College 

We  are  studying  methods  of  reconstructing  the  elasticity  from  MRI  measurements 
of  tissue  vibration.  There  has  been  significant  interest  in  estimating  tissue  elasticity  from 
MRI  phase  contrast  measurements  of  periodic  and  quasi-static  displacement.  MRI  seems 
to  hold  more  promise  than  ultrasound  because  of  its  ability  to  measure  small  tissue 
displacements  simultaneously  in  all  three  directions  resulting  from  a  single  mechanical 
stimulus  while  ultrasound  is  limited  to  recording  tissue  displacements  in  one  preferred 
direction  at  at  time. 

We  have  calculated  tissue  displacements  with  the  partial  differential  equations 
describing  dynamic  and  static  elastic  deformation.  Models  of  the  breast  were  generated 
from  MRI  scan  data.  We  have  performed  simulations  for  various  modes  of  vibration. 
These  simulations  have  led  to  three  conclusions  which  impact  how  estimates  of  elasticity 
can  be  obtained  from  displacement  fields: 

1)  If  the  driving  displacement  is  large  enough  to  obtain  3D  MR  phase  contrast  images  in 
reasonable  times,  there  is  likely  to  be  significant  displacement  in  directions  perpendicular 
to  the  direction  of  the  driving  force. 

2)  Multi-dimensional  displacement  (e.g.  in  directions  other  than  in-line  with  the  driving 
force)  requires  partial  differential  equation  solution  to  adequately  describe  the 
displacement  field. 

3)  Because  partial  differential  equations  are  necessary  to  describe  the  motion,  those 
equations  must  be  used  to  estimate  the  elasticity. 

If  the  displacement  is  all  essentially  in  the  direction  of  the  driving  force,  simple  local 
estimates  of  the  elasticity  would  be  possible. 
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Introduction:  ■  ■  c 

We  introduce  a  new  method  of  removing  noise  trom 
images.  The  method  selects  significant  maxima  (includes 
the  minima)  in  the  image  and  forces  the  pixel  values  in  the 
filtered  image  to  vary  monotonically  in  all  directions 
between  those  maxima.  For  example,  if  an  isolated 
maximum  is  identified,  the  filtered  image  should  fall  away 
from  it  monotonically  in  all  directions.  We  can  only 
perform  monotonic  fits  in  one  dimension,  so  we 
approximate  monotonicity  in  all  directions  by  doing 
monotonic  fits  along  line  segments  through  out  the  image. 
The  line  segments  are  bounded  by  maxima  and  are  onented 
in  many  directions  through  the  image.  The  filtering 
operation  on  each  line  segment  replaces  the  pixel  va^es  on 
that  segment  with  a  ID  monotonic  sequence  that  fits  the 
original  pixel  values  best  in  a  least  squares  sense.  The 
method  is  simple,  relatively  fast  and  stable.  SNR’s  as  low 
as  0.1  can  produce  acceptable  results. 

Methods:  .  ^  j  u  i 

An  intermediate  filtered  image  is  found  at  each  angle. 
The  original  image  is  rotated  to  the  appropriate  angle.  The 
rotated  image  is  blurred  using  a  simple  Gaussian  filter  in 
the  Fourier  domain.  The  maxima  along  the  rows  of  the 
image  are  identified.  The  maxima  that  differ  from  the 
adjacent  maxima  by  less  than  the  threshold  are  removed 
from  the  set.  In  each  row,  all  sequences  of  adjacent  pixels 
ending  in  maxima  are  replaced  by  the  monotonic  sequence 
that  best  fits  the  original  pixel  values  [1,2].  Then  the  same 
procedure  is  followed  on  the  columns  of  the  row  filtered 
image.  Remarkably,  the  rows  stay  essentially  monotonic 
when  the  columns  are  processed.  The  row-column  filtered 
image  is  averaged  with  the  column-row  filtered  image  and 
rotated  back  to  the  original  orientation.  The  result  is  the 
intermediate  filtered  image  for  that  angle.  Int^ediate 
filtered  images  are  obtained  for  each  angle.  The  final 


filtered  image  is  a  least  squares  combination  of  all  of  the 
intermediate  filtered  images.  The  final  filtered  image  is 
almost  monotonic  between  maxima  in  all  directions. 
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Fig.  1  An  example  of  the  filter  on  a  row.  a)  signal:  two 
boxcars  with  noise,  SNR«1.  b)  the  blurred  signal  and  the 
identified  maxima,  c)  result  of  the  monotonic  fits 
superimposed  on  the  original  boxcars.  Peaks  in  the 
filtered  signal  is  noise  that  passes  through  the  filter  at  the 
maxima  when  the  SNR  is  very  low.  Note  that  the  edges  of 
the  recovered  boxcars  are  sharp  &  in  the  correct  positions. 

ronclusions:  ^ 

We  have  presented  a  method  of  removing  noise  from 
images  using  monotonic  fits  between  maxima.  It  does  not 
blur  edges,  is  robust  and  can  be  effective  on  images  with 
very  low  SNR’s  unlike  most  wavelet  denoising  methods. 
However,  it  does  not  yet  have  the  mathematical  structure 
and  rigor  that  wavelet  denoising  methods  have.  It  also 
leaves  point  noise  at  the  maxima  at  very  low  SNR  s. 

1  I.C.  Demetriou  and  M.J.D.  Powell,  IMA  Journal  of  Numerical 

Analysis \\'A\\-Ayi  \99\.  oTin -70  loo^ 

2  J.B.  Weaver,  et  al.  SPIE:  Medical  Imaging  2710-79,1996. 


Fie  2  a>  orifiinal  high  SNR  image,  SNR=80  within  brain,  b)  image  with  added  noise,  SNR- 16.  c) 

®d)  gLS  bluSl  b  sitowini  (he  maxima,  e)  result  of  filter  on  b,  SNR=85.  0  result  of  filter  on  c,  SNR=30. 
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Abstract 

We  have  developed  a  method  of  reconstructing  the  tissue  e- 
lasticity  from  MR  displacement  measurements.  The  model 
used  to  reconstruct  the  elasticity  is  based  on  the  partial  dif¬ 
ferential  equations  describing  steady  state,  harmonic  motion. 
Therefore,  we  have  developed  a  system  to  generate  true  steady 
state  motion  within  tissue  as  well  as  within  gel  samples.  Our 
first  motion  maps  and  reconstructions  from  a  homogeneous  gel 
phantom  are  shown. 

Introduction 

Measurements  of  harmonic  motion  and  quasi-static  displace¬ 
ment  can  be  accomplished  accurately  using  MR  [1,2,3].  The 
inversion  of  that  data  to  estimate  the  mechanical  properties 
of  the  tissue  has  proved  to  be  a  difficult  problem.  Currently 
several  groups  have  taken  different  approaches  and  the  MR 
acquisition  has  been  modified  to  fit  the  inversion  approach. 
Methods  to  estimate  the  elasticity  from  traveling  waves  and 
from  quasi-static  displacements  have  been  developed.  In  all  of 
these  methods,  reconstruction  has  been  problematic.  We  have 
implemented  an  inversion  method  using  finite  element  solu¬ 
tions  to  the  partial  differential  equations  describing  the  elasto- 
dynamic  response  of  soft  tissue  under  an  applied  harmonic  de¬ 
formation.  The  inversion  method  is  a  nonlinear  Newt  on- based 
iterative  scheme  which  finds  the  distribution  of  Young’s  modu¬ 
lus  which  minimizes  the  difference  between  the  measured  and 
model-predicted  displacement  values  [4].  The  problem  with 
our  previous  MR  acquisitions  is  that  the  motion  is  started  and 
stopped  during  each  cycle  of  the  pulse  sequence  and  we  have 
suggested  that  the  transients  in  the  motion  produce  incorrect 
reconstructions  [5,6].  Therefore,  we  have  developed  a  method 
of  measuring  displacements  during  steady  state  motion. 

Methods 

With  simulations  we  have  estimated  that  100  cycles  of  100 
Hz  motion  allow  the  transients  to  dissipate  for  a  wide  variety 
of  mechanical  properties  and  geometries.  It  is  impractical  to 
get  100  cycles  of  motion  for  each  MR  excitation  so  we  kept 
the  motion  going  throughout  the  gradient  echo  pulse  sequence. 
The  motion  was  small  enough  that  it  did  not  disturb  the  signal 
significantly  except  when  the  motion  sensitizing  gradients  were 
on.  The  difference  between  an  image  with  and  without  the 
motion  was  negligibly  small  without  the  motion  sensitizing 
gradients,  which  were  two  sine  waves  just  before  the  read-out. 

The  signal  used  to  drive  the  gel  sample  was  produced  by  an 
HP  331 20 A  signal  generator  that  was  phase  locked  to  the  10 
MHz  MR  system  clock.  The  signal  generator  was  setup  to  gen¬ 
erate  a  set  number  of  100  Hz  sine  waves  after  it  was  triggered 
with  a  pulse  from  the  MR.  The  signal  generator  then  drives  the 
sample  at  100  Hz  for  the  duration  of  the  pulse  sequence.  The 
initial  phase  of  the  100  Hz  signal  was  controlled  by  the  signal 
generator  and  could  be  set  by  the  operator.  That  initial  phase 
controls  the  relative  phase  between  the  motion  and  the  motion 
sensitizing  gradients  which  must  be  changed  to  determine  the 
motion  in  each  voxel.  To  keep  the  motion  synchronized  with 
the  pulse  sequence  throughout,  the  actual  TR  must  be  an  in¬ 
teger  multiple  of  the  10  msec  cycle.  The  measured  TR  was 
actually  1  msec  more  than  the  TR  set  on  the  console.  The 
signal  from  the  signal  generator  was  the  input  of  a  power  am¬ 


plifier  that  drove  a  set  of  6  piezoelectric  stacks  that  vibrated 
the  sample. 

Results 

Our  first  maps  of  steady  state  motion  in  a  homogeneous  gel 
sample  are  shown  below.  The  reconstruction  of  the  elasticity 
is  also  presented.  The  average  elasticity  in  the  image  is  ap¬ 
proximately  what  we  expect  from  physical  measurements  in 
the  gel,  around  20  kPa.  The  reconstruction  is  fairly  unifor- 
m,  although  there  is  some  structure  in  the  image  where  the 
recovered  values  deviate  from  the  nominal  level. 


X  Phase  Y  Phase 


Figure  1:  The  amplitudes  and  phases  of  the  harmonic  steady 
state  motion  in  an  essentially  homogeneous  agar  gel  phantom. 


Figure  2:  The  elasticity  reconstructed  from  the  motion  shown 
above  in  an  essentially  homogeneous  agar  gel  phantom. 

This  work  was  supported  by  NIH  grant  R01-NS33900,  NSF 
grant  BCS-9978116  and  NIH  grant  POl  CA80139-01A1. 
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Abstract 

Recently,  a  variety  of  methods  for  generating  elastic  con¬ 
trast  images  of  biological  tissue  have  been  put  forward  [1,2,3]. 
To  date  the  large  majority  of  these  methods  have  been  devel¬ 
oped  under  assumptions  of  two  dimensional  mechanical  behav¬ 
ior  for  the  tissue  in  question.  While  this  may  well  be  valid  in 
certain  idealized  conditions,  it  is  by  no  means  accurate  for  the 
general  case,  especially  for  non-symmetric  tissue  geometries 
or  property  distributions.  In  these  situations,  a  three  dimen¬ 
sional  property  reconstruction  scheme  will  be  necessary.  Here, 
our  efforts  to  develop  a  three  dimensional  elastic  property  re¬ 
construction  scheme  are  discussed  and  simulation  results  are 
presented. 

Introduction 

MR  based  elastic  property  imaging  is  a  rapidly  evolving  field 
which  seeks  to  determine  the  mechanical  property  distribution 
within  a  tissue  region  using  displacement  or  strain  informa¬ 
tion  obtained  for  that  region.  This  approach  is  necessary  as 
a  tissue’s  characteristic  elastic  makeup  is  not  directly  visible 
by  means  of  magnetic  resonance  or  ultrasound  imaging.  For 
the  most  part,  the  schemes  developed  for  reconstructive  imag¬ 
ing  have  been  presented  using  two  dimensional  assumptions. 
While  such  assumptions  are  valid  in  certain  cases,  namely  per¬ 
fect  symmetry  in  both  property  distribution  and  geometry, 
they  will  not  be  valid  in  a  general  heterogeneous  tissue,  such 
as  the  human  breast.  In  these  instances  motion  or  strain  in 
the  third  dimension  is  not  insignificant  and  must  be  accounted 
for. 

Methods 

Previously,  a  finite  element  implementation  of  a  subzone 
based  elasticity  reconstruction  scheme  has  been  presented  [2]. 
This  method  uses  the  full  field  displacement  data  available 
from  the  MR  to  drive  a  non-hnear  reconstruction  process  based 
on  squared  error  minimization.  This  approach,  documented 
using  a  two  dimensional  plane  strain  approximation,  is  fully 
amenable  to  a  complete,  three  dimensional  treatment.  Using  a 
linear  elastic  model  as  the  basic  assumption  for  tissue  motion, 
the  governing  equation  for  the  harmonic  tissue  response  is 


V  •  GVu  +  V(A  +  G)V  -  u  =  pJ-a.  (1) 

The  inversion  process  based  upon  the  three  dimensional 
equations  of  linear  elasticity  is  executed  in  a  similar  fashion 
to  the  two  dimensional  scheme,  with  one  notable  exception. 
The  automated  zoning  process  for  the  two  dimensional  inver¬ 
sion  problem  is  driven  by  a  hierarchical  ordering  of  element 
based  error,  requiring  a  global  solution  using  the  current  prop¬ 
erty  distribution  estimate  for  every  global  iteration.  For  high 
resolution  three  dimensional  problems  global  solutions  are  too 
costly  to  perform  at  each  global  inversion  iteration.  Thus,  for 
the  three  dimensional  inversion  algorithm  the  zone  location  is 
selected  in  a  random  manner  from  the  list  of  elements  not  yet 
operated  on  during  the  current  iteration  step. 

Results  k,  Discussion 


using  displacements  generated  through  a  finite  element  solu¬ 
tion  with  15%  random  noise  added  to  the  data.  To  minimize 
the  effects  of  such  added  noise  on  the  inversion  solution,  a 
small  degree  of  spatial  filtering  was  used  during  the  inversion 
process.  This  is  achieved  by  incorporating  the  current  prop¬ 
erty  estimate  of  nodes  in  direct  contact  with  a  given  node,  z, 
such  that  Er'“  =  (1  -  e)Ef'^  +  where  N,  is 

the  number  of  neighbor  nodes  connected  to  node  i.  For  this 
experiment  a  value  of  0.05  was  used  for  0. 


Figure  1:  Three  dimensional  phantom  simulation  reconstruc¬ 
tion  for  a  5cm  x  5cm  x  3.5cm  geometry  with  a  25  kPa  back¬ 
ground  Young’s  modulus  and  a  1  cm  diameter  250  kPa  spher¬ 
ical  inclusion,  (a)  Inversion  solution  for  0%  noise  case,  (b) 
Inversion  solution  for  15%  noise  solution. 


The  inversion  process  consisted  of  roughly  20  global  itera¬ 
tions,  each  consisting  of  300  zone  based  parameter  updates  on 
average,  with  the  average  zone  size  being  roughly  330  nodes. 
The  total  solution  consists  of  24094  nodal  stiffness  values.  This 
work  was  supported  in  part  by  NIH  grant  530374. 
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Initial  simulation  experiments  have  been  performed,  fig.  1, 


